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Abstract 


o/An  effective  numerical  method  for  solving  boundary  value  problems  for  the  Landau 
Fokker -Planck  equation  is  developed  and  applied  to  calculating  the  electron  velocity 
distribution  function  in  model  solar  transition  regions.  Numerical  results 
illustrating  the  speed,  pitch  angle  and  spatial  dependence  of  the  distribution 
function  are  presented.  From  these  it  is  concluded  that  the  widely-invoked  assump¬ 
tion  that  in  weakly  inhomogeneous  colli sional  plasma  the  angle -averaged  distribution 
function  remains  close  to  local-Maxvellian  distribution  is  Incorrect.  Instead, 
the  distribution  function  forms  an  anisotropic,  high  velocity  tail  In  the  lower 
temperature  regions  due  to  the  diffusion  of  fast  electrons  anti-parallel  to  the 
temperature  gradient .  I  Roughly  speaking  this  effect  is  of  quantitative  slgnflcance 


for  spectroscopic  work  and  heat  conduction  provided  1  *  10”J,  where 

2  4  y 

X  •  (kT)  /ire  nipAr  is  an  effective  mean  free  path  for  thermal  electrons. 

^Et  is  shown  that  as  a  result  of  there  being  an  excess  of  fast  electrons  in  the 

5^  peiver' 

low  transition  region  (T  *  3  x  K,  say) ,  inelastic  electron-ion  collision  rates 

A  A 

are  significantly  enhanced  over  the  Maxwellian  values.  Attendant  effects  on  the 
ionization  balance  of  a  typical  metal  (magnesium)  are  shown  to  be  significant. 
Implications  of  the  breakdown  of  the  local-Maxwelllan  approximation  for  several 
outstanding  questions  related  to  the  solar  transition  region  are  discussed.  Including: 
energy  balance  In  low  transition  region  and  upper  chromosphere,  the  helium  resonance 
line  spectrum,  the  Sehmahl-Orrall  observation  of  continuum  absorption  by  neutral 
hydrogen,  and  the  origin  of  the  20,000  K  temperature  plateau. A 
Subject  headings:  plasmas  -  Sun:  atmosphere  -  ! 

i 

Sun:  chromosphere  -  atomic  processes.  ' 


i  10~^,  where 


I .  Introduction 

This  is  ths  first  in  a  projected  series  of  papers  on  the  topic  of  how  gradien 
in  density  and  temperature  affect  the  electron  velocity  distribution  function 
(EVDF)  in  collislonal  plasma.  For  reasons  to  be  discussed,  emphasis  through¬ 
out  the  series  will  be  placed  on  physical  conditions  thought  to  exist  in  the 
Sun's  chromosphere-corona  transition  region  and  upper  chromosphere  (hereafter 
referred  to  simply  as  the  transition  region;  TR)  and  on  assessing  any 
implications  of  the  kinetic  theory  to  be  developed  for  the  spectroscopic 
diagnosis  and  energy  balance  of  these  regions.  Results  of  more  general 
Interest  will  also  be  presented,  however. 

Knowledge  of  the  form  assumed  by  the  EVDF  in  a  plasma  under  weakly 
inhomogeneous  conditions  -  those  in  which  the  effective  mean  free  path  of 
thermal  electrons  X  (X  •  (kT)a  /xe*nlnA  in  fully  ionized  hydrogen;  see  below) 
is  everywhere  small  relative  to  gradient  scale  lengths  -  is  an  essential 
element  in  many  astrophysical  analyses.  Such  knowledge  is  required,  for 
example,  in  order  to  determine  electron-ion  collislonal  excitation  and 
ionization  rates  and  the  plasma  thermal  and  electrical  conductivities. 

Although  the  Sun's  TR  is  known  from  EUV  observations  to  be  the  site  of 
unusually  large  temperature  and  density  gradients,  it  nonetheless  qualifies 
as  a  weakly  inhomogeneous  medium  under  the  definition  given  above.  This 
point  is  illustrated  by  Table  1  showing  portion  of  an  empirically  derived, 
spatially  averaged,  quiet-sun  TR  model  published  by  Dupree  (1972.)  Note 
that  a  =  X  <  5  x  10_*.  Consequently,  it  is  generally  thought  that 

throughout  the  Sun's  TR,  as  in  any  quasi-stationary ,  weakly  inhomogeneous 
plasma,  the  EVDF  can  be  adequately  approximated  by  an  expression  of  the 


(i-,  v,  z)  -  f*  (v;  n,  T)  j  l-H:  i(z)l>  (-2—  \  +  O(ft’)  ! 


where  I  have  considered  the  simplest  case  of  a  horizontally  stratified 

(VT  *  —  e  ),  isobarlc  plasma  at  rest  with  a  vertical  (or  zero)  magnetic 
~  dz  z 

field.  Here  f *  Is  a  local-Maxwellian  distribution  with  number  density  n 

v.4 

and  temperature  T,  y  •  ,  and  D  (x)  is  obtained  by  applying  the 

Chapman-Enskog  formalism  to  a  particular  kinetic  equation,  as  done,  for 

example,  by  Spitzer  and  Harm  (1953)  for  the  Fokker-Planck  equation. 

Equation  (1)  has  two  important  implications:  first  the  angle^averagad 

distribution  f  *  ^  /*  fdy,  remains  equal  to  a  local  Maxwellian  distribution, 

00 

and  second,  the  electron  heat  flux,  q  =  2ff^./lydy  f  v5  f  (y,  v,  z)  dv,  is 
5 /2  dT  ^ 

proportional  to  T  —  and  is  only  logarithmically  dependent  on  density.  These 
conclusions  underlie  many  solar  and  astrophysical  analyses,  including  ionization 

equilibrium  and  radiative  loss  calculations,  spectroscopic  diagnostic 
procedures  and  energy  balance  arguments.  It  is  my  aim  in  this  and  later 
papers  in  this  series  to  show  that  they  are,  in  fact,  invalid  in  the  solar 
transition  region  -  and,  by  implication,  under  similar  conditions  in  other 
inhomogeneous  media. 

The  difficulty  is  that  equation  (1)  is  not  a  uniformly  valid  approxima¬ 
tion  (Shoub  1976,  1981);  it  fails  for  electron  velocities  v  >  v  ,  where  v 

/vc  C  C 

is  defined  bv  the  relation  a  ( - 1  ■  1.  The  reason  for  this  failure  is 

\Vth/ 

that  according  to  Landau* s  (1936;  1949)  Fokker-Planck  equation  (the  same 
equation  was  later  derived  using  different  arguments  by  Spitzer 
and  Harm  (1953)  and  later  still  by  Rosenbluth, 


MacDonald  and  Judd;  RMJ  1957),  the  effective  speed -dependent  mean  free  path 
of  an  electron  in  a  fully  ionized  gas  increases  as  the  fourth  power  of  its 
velocity,  for  velocities  greater  than  thermal.  As  defined  above,  v  is  the 


velocity  at  which  an  electron's  mfp  equals  the  local  temperature-gradient 


scale  length.  The  ratio  v  /v  .  is  listed  in  Table  1;  note  that  it  falls 

c  tn 

well-within  the  range  relevant  to  calculating  inelastic  collision  rates. 
Equation  (1)  fails  for  v>v^  because  its  derivation  is  based  on  the  assumed 
existence  of  "normal",  or,  in  present  context,  spatially  local  solutions  to 
a  kinetic  equation  of  the  form 

f  (v,  z)  «  f  (v;  n(z) ,  T(z) ,  Vn(z),  ?T,  ...) 
where  the  ellipsis  denotes  higher  order  gradients.  Zt  seems  clear,  at  least 
intuitively,  that  such  a  local  solution  cannot  exist  at  velocities  for  which 
electrons  can  traverse  one  or  more  scale  heights  before  thermalizing.  This 
can  be  seen  from  a  mathematical  standpoint  as  follows. 

Consider  the  dimensionless,  high-velocity  form  of  the  RMJ  Fokker-Planck 
equation  in  the  Haxwellian-field-partlele  approximation  (hereafter  MFP  A) , 


namely 

vC 


+ 

X  + 


jL  !(i-u2y& 

3u  U  V  oU 


d*n(TQ) 

dr 


Equation  (2)  follows  from  the  full  Fokker-Planck  equation  upon  evaluating  the 
Rosenbluth  potentials  (see  Eqn.  9)  using  a  local-Maxwellian  distribution  for 
electrons,  a  delta  function  distribution  at  zero  velocity  for  protons  and 
then  neglecting  terms  in  the  resulting  coefficients  which  are  unimportant 


at  large  velocities.  A  thermal  electric  field  term  has  also  been  omitted 

from  the  left  side  of  (2);  it  is  unimportant  relative  to  gradient-related 

terms  at  large  velocities  (see  Eqn.  41).  The  quantities  T  and  n  are 

o  o 

reference  variables  which  for  the  moment  may  be  identified  with  the 
electron  temperature  and  density..  Note  that  X,  referred  to  above  as  the 
effective  mean  free  path  of  thermal  electrons,  is  simply  the  length  scale 
which  emerges  from  non-dimens localizing  the  Fokker-Planck  equation.  Its 

numerical  value  depends  on  the  choice  of  characteristic  velocity  v 

(v  \* 

—  |  .  The  left  side  of  Eqn.  (2)  is  proportional 

vth/ 

to  and  incorporates  the  relation  n  T  -  constant.  The  5 -derivative 
az  oo 

terms  on  the  right  derive  from  electron-electron  collisions  while  the 

y-derlvative  term  contains  equal  contributions  from  electron-electron  and 

electron-proton  scattering.  A  factor  of  £*  has  been  brought  to  the  left 

side  in  (2)  to  facilitate  the  following  argument.  Note,  finally,  that 

2  -£2 

disregarding  boundary  conditions,  a  local-Maxwellian  <J>*  ■  ,  satisfies 

Eqn.  (2)  when  oc(t)  =  0.  Eqn.  (2)  is  derived  in  detail  in  SlI. 

Consider  now  the  traditional  (Chapman-Enskog)  argument  leading  from 

Eqn.  (2)  to  Eqn.  (1).  First,  one  seeks  a  spatlally-local  (normal)  solution 

of  (2)  in  which  there  is  no  explicit  spatial  variation  of  but  only  an 

implicit  one  arising  from  its  dependence  on  the  thermodynamic  variables  and 

their  gradients.  In  the  context  of  Eqn.  (2)  this  means  we  look  for  solutions 

3a 

of  the  form  <j>  (u,  5,  t;  a(x))  ■  <Ku,  £»  *»  »  •  .  •).  Second,  since 

a  is  presumed  small,  one  assumes  that  <£  nay  be  expanded  as  a  power  series  in 
the  a.  (This  is  known  in  the  kinetic  theory  literature  as 
the  Hilbert  expansion  assumption;  see,  e.g.,  Resibos  and  De  Leener  1977). 


$  :  v°(y,  5)  +  £)  +  .  . 


On  substituting  this  expansion  into  (2)  and  noting  that  *  0  (a2) , 

one  finds,  in  the  usual  way,  that 

$(W,  5,  T)  -  <j>*(CHl  +  ya  $*($)  +  o  (a2)}, 
where  is  a  local -Maxwellian  and  41  satisfies  an  inhomogeneous  ordinary 
differential  equation.  This  procedure,  then,  is ,  in  essence,  the  origin  of 
the  well-known  Spitzer-Harm  (1953)  solution.  It  is  seen  that  the  approxima¬ 
tion  of  local  thermodynamic  equilibrium  (LTE)  -  taken  here  to  mean  that 
—  f  1 

2J  4  dli  »  ♦*  -  is  a  necessary  consequence  of  the  assumption  that  the 
distribution  function  may  be  expanded  in  a  power  series  in  the  gradients, 
for  a  local-Maxwellian  must  be  annihilated  by  any  viable  collision  integral. 
LTE  is  therefore  built  into,  rather  than  deduced  from,  the  Chapman-Enskog 
analysis. 

But  inspection  of  Eqn.  (2)  clearly  shows  that  the  regular-perturba¬ 
tive  approach  just  described  cannot  give  valid  results  at  velocities  for 
which  a£4  j>  0  (1).  In  this  regime  the  gradient-related  terms  in  (2)  dominate 
the  collision  terms,  and  this  precludes  finding  meaningful  spatially-local 
solutions. 

The  actual  form  assumed  by  the  tail  of  the  EVDF  under  conditions  for 
which  a£4  becomes  of  order  unity  at  moderate  velocities  is  thus  an  open 
question  -  and  is  the  question  addressed  in  this  paper.  The  approach  taken 
has  been  to  mount  a  careful  numerical  attack  on  the  RMJ  Fokker-Planck 
equation  -  in  the  present  context  a  quasl-linear,  second-order,  partial 
differential  equation  in  three  independent  variables  with  half-range  boundary 
conditions.  The  numerical  method  developed  is  a  fourth-order,  fully  implicit 
finite  difference  algorithm.  After  extensive  testing,  the  code  was  applied 


Co  calculating  Che  EVDF  for  Che  simplest  relevant  models  of  the  solar  TR, 
namely  inhomogeneous  plasma  slabs  of  fully  ionized  hydrogen  with  a  vertical 
(or  zero)  magnetic  field.  (Optically  thin  radiation  is  also  allowed  for, 
but  not  self-consistently.)  Inhomogenlty  is  introduced  via  boundary 
conditions  on  the  incoming  distribution  functions.  The  main  conclusion 
drawn  from  these  calculations  is  that  the  local-Maxwellian  approximation  is 
badly  in  error  throughout  the  middle  and  lower  TR.  The  computed  distributions 
exhibit  a  pronounced,  anisotropic,  high-velocity  tail  attributable  to  the  free 
streaming  of  fast  electrons  from  hotter  overlying  layers.  This  effect  causes 
collislonal  ionization  rates  (and  in  some  cases  excitation  rates)  to  be 
enhanced  over  their  LTE-values,  often  by  several  orders  of  magnitude. 

The  numerical  results  reported  here  have  been  obtained  in  the  Maxwell ian- 
f  leld-par tide-approximation  (MFPA.) ,  meaning  that  the  Rosenbluth  potentials 
which  enter  the  Fokker-Planck  equation  are  evaluated  using  a  local-Maxwellian 
distribution  for  electrons.  In  IIII  I  show  that  although  this  approximation 
gives  the  angle-averaged  part  of  the  distribution  function  accurately,  its 
use,  together  with  my  not  having  calculated  the  distribution  function  out  to 
sufficiently  large  velocities,  precludes  an  accurate  determination  of  the 
electron  flux.  The  heat  flux  question  is  therefore  deferred  to  a  later 
publication. 

The  rest  of  the  paper  is  organized  as  follows:  In  SlI  I  present 
the  mathematical  model  and  discuss  its  limitations.  In  §111  I  discuss  the 
numerical  algorithm  and,  in  §IV,  give  detailed  results  for  the  distribution 
function,  showing  its  spatial,  velocity  and  pitch  angle  dependence.  In  this 
section  I  also  calculate  inelastic  collisi  ;n  rates  for  helium  and  magnesium 
and  show  how  the  enhanced  ionization  rates  alter  magnesium 41  s  ionization 
equilibrium.  In  §V  I  discuss  possible  Implications  of  the  results  and 


suggest  chat  several  observational  and  theoretical  puzzles  concerning  the 
low  TR  and  upper  chromosphere  might  be  resolved  by  consideration  of  this 


W" 


8 


F* 


II .  Model  Definition 

I  consider  an  Idealized  TR  consisting  In  a  constant  pressure  slab  of 

fully  ionized  hydrogen  containing  a  fictitious  trace  ion  in  sufficient 

amount  to  allow  the  plasma  to  radiate  energy  at  a  rate  equal  to  that 

radiated  by  an  optically  thin  plasma  with  cosmic  element  abundances,  as 

calculated  by  McWhirter  £t  al.  (1975) .  The  slab  has  thickness  L  and  any 

magnetic  field  present  is  assumed  to  lie  along  the  gradient  direction.  The 

protons  are  taken  to  be  infinitely  massive,  at  rest,  and  distributed  so  as  to 

provide  charge  neutrality.  I  neglect  the  trace  ion's  contribution  to  the 

charge  balance  so  that  n  ■  n  =  n.  The  constraint  that  no  steady  current 

P  e 

flow  through  the  TR  then  implies  that  the  average  electron  velocity  must  be 
zero.  Inhomogenity  is  introduced  via  boundary  conditions  on  the  Incoming 
electron  distribution  function  on  the  planes  z  *  0  and  z  »  L,  as  described 
below.  Finally,  I  assume  that  the  EVDF  obeys  the  Fokker-Planck  equation, 
in  the  form  derived  by  Rosenbluth,  MacDonald  and  Judd  (1957) .  Several 
limitations  of  this  central  assumption  should  be  mentioned. 

First,  the  RMJ  equation  is,  at  best,  accurate  only  to  within  terms  of 
order  (in  A)  1 ;  i.e.,  to  "dominant  order"  (e.g.,  Grad  1962).  Moreover, 
it  appears  (Spitzer  1962,  Siambis  and  Stitzer^  1974)  that  the  "dominant" 
approximation  is  not  uniformly  valid  in  velocity.  More  specifically,  a 
non-dominant  contribution  to  the  parallel ' (ev  §v)  component  of  the  velocity 
diffusion  tensor,  which  of  course  is  omitted  from  the  RMJ  equation,  decays 


1  Siambis  and  Stitzer  erroneously  attribute  differences  between  their 

results  for  the  parallel  diffusion  coefficient  and  those  of  Spitzer  (1962) 
(and  RMJ)  to  their  retention  of  a  velocity-dependent  Coulomb  logarithm. 

The  latter  refinement  leads  only  to  minor  corrections  to  the  Spitzer 
result  however.  The  main  reason  for  the  discrepancy  lies  in  their 
keeping  a  non-dominant  term  which  is  neglected  by  Spitzer. 


9 


as  v  1  for  large  v,  whereas  Che  dominant  contribution  decays  as  v”* .  The 

neglected  term  is  therefore  no  longer  negligible  for  )  2  >0  (in  A). 

*  vth  / 

Nevertheless,  Siambis  and  Stitzer  show  that  the  contribution  of  the  neglected 
term  to  the  energy  loss  rate  of  a  test  particle  is  of  relative  order  (in  A)-1 
at  all  velocities,  thereby  leading  one  to  suspect  that  its  retention  would 
not  significantly  alter  results  obtained  from  the  RMJ  equation.  This  latter 
conjecture  is  currently  under  investigation. 

Second,  because  the  RMJ  equation  was  derived  from  consideration  of 
successive  but  uncorrelated  binary  events,  its  use  precludes  consideration 
of  possible  collective  effects.  However,  in  the  present  problem,  where  the 
question  is  whether  significant  departures  from  a  local-Maxwellian  distribu¬ 
tion  exist  under  conditions  in  which  classical  theory  suggests  they  do  not, 
the  disregard  of  collective  effects  may  be  viewed  as  a  plausible  working 
hypothesis,  which  may  be  checked  a  posteriori,  at  least  in  principle.  In 
practice,  furnishing  proof  that  collective  effects  are  unimportant  is  likely 
to  be  difficult.  Not  only  must  one  demonstrate  stability  against  both 
electrostatic  and  electromagnetic  perturbations,  but  one  must  also  show  that 
electron  scattering  by  enhanced  plasma  fluctuations,  which  arise  even  in 
stable  non-Maxwe Ilian  plasmas  (e.g.,  Tidman  and  Eviatar  1965),  is  negligible 
compared  to  Coulomb  scattering.  Clearly,  then,  the  question  of  collective 
effects  in  the  present  context  is  an  involved  one,  and  is  likely  to  be  a 
fruitful  area  of  future  research.  It  is  not  addressed  further  in  this  paper. 
The  reader  should  bear  in  mind  that  such  effects  may  alter  conclusions 
drawn  below. 

Third,  and  lastly,  note  that  because  the  effect  of  a  background 
magnetic  field  on  binary  collision  dynamics  is  ignored  in  the  RMJ  equation, 
its  use  is  limited  to  situations  in  which  the  mean  electron  gyroradius 


10 


is  large  compared  with  the  Debye  length  (e.g. ,  Baldwin  and  Watson  1975).  This 
-2  1/2 

requires  B  <  10  n  ,  where  B  is  the  magnetic  field  strength  in  gauss 

-3 

and  n  is  the  electron  density  in  cm  .  This  condition  is  not  overly 

restrictive  for  solar  applications. 

Consider  now  the  appropriateness  of  a  time-independent  calculation. 

-14  -3 

For  a  TR  with  nT  ■  const.  *  6  x  10  K  cm  ,  the  electron-electron  relaxa- 
cion  time,  v_1(v)  -  —  (—)  3  y  3.5  *  1<T16 

“  3  Vch  Khj  ' 

5  x  10  ^  (—^—\  sec  at  the  top  of  the  TR  (T  ■  10^*®  K,  in  A  ■  20)  and 


(vch) 


T  <—\ 
A  rth  I 


,  varies  between 


3  x  10  ^  V  sec  at  the  bottom  (T  »  10^  K,  in  A  ■  11).  Thus,  even  in 


te”) 


the  far  tail  of  the  distribution  (  ^51,  stationarity  is  achieved  in 

_4 

approximately  one  second  at  the  top  of  the  TR  and  in  10  sec  at  the  bottom. 
On  the  other  hand,  relevant  macroscopic  time  scales,  for  example,  the  rise 
times  of  spicules,  are  typically  of  the  order  of  minutes  (e.g.,  Athay  1976,  p. 
Proceeding,  I  write  the  electron  Fokker-Planck  equation  as 


af  e  v  v  af  6f 

v  •  f  -  f  <E  ♦  ?  *  B)  •  €  ^ 


el 


6f 

ot 


inel. 


(3) 


Here  -e  and  m  are  the  electron  charge  and  mass,  E  is  the  self-consistent 
polarization  electric  field  required  for  zero  current  flow,  and  the  terms 
on  the  right  represent  elastic  and  inelastic  collision  terms.  The  latter 
are  included  to  allow  for  the  effect  of  radiative  losses  on  the  electron 
temperature.  The  distribution  function  f  is  normalized  to  the  average 
electron  number  density. 

Let  (e  ,  e  ,  £  )  be  an  orthonormal  triplet  with  §  in  the  direction  of 
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increasing  temperature.  I  adopt  a  spherical  coordinate  system  in  velocity 
space  (v,  0,  8)  with  polar  axis  along  8^.  The  corresponding  unit  vectors  are 
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cos  ?  e  +  sin  8  (cos  8  e  +  sin  Be) 
z  x  y 


-  sin  8  £  +  cos  8  (cos  M  +  sin  Be) 

Z  V  v 


■  -  sin  M  +  cos  B  £ 

B  x  y 


It  will  be  useful  below  to  note  that 


n  A  3  «  1  3  8  3 

7v  “  \  57  +  *8  v  38  +  v" sin  0  36 


and  that 
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with  all  other  unit-vector  derivatives  equal  to  zero. 


Now  consider  the  magnetic  field  terms  In  (3).  With 


Be  I  find 
z 


e  ,  _ .  .  3f  n3f 

5^*n38’ 


where  C  •  eB/mc  is  the  electron  gyrofrequency .  The  spatial  gradient  term 


in  (3)  is 
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v  •  t—  =  v  cos  8  t—  +  sin  9 
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where  the  second  equality  follows  from  the  assumed  homogeneity  of  f  in 
the  x-y  plane.  If,  in  the  boundary  conditions, 

f(v,  2  -  0)  -  f£  (y)  J  v  *  %z  >  0 

(8) 

f(v,  2  -  L)  -  f“  (v)  ;  V  •  ez  <  0, 

the  prescribed  functions  f^  are  chosen  8-independent ,  it  then  follows  that 
3f 

g-g  =  0,  for  there  are  no  further  terms  in  (3)  which  introduce  a  B -dependence 
into  f .  Thus,  under  the  above  restrictions  the  magnetic  field  has  no 
effect  on  the  distribution  function.  'Viewed  differently,  a  vertical  magnetic 
field  is  irrelevant  because,  for  a  given  pitch  angle,  the  helical  and 
slant  pathlength  between  heights  s^  and  are  equal. 

Following  Roseribluth  et  al.  (1957)  I  write 


(9a) 


(9b) 


ha(z»v) 


ga(z,v) 
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(9c) 


(9d) 
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Here  ur.subscrlpted  species-dependent  quantities  refer  to  electrons  and 
the  summation  in  (9a)  is  over  electrons  and  protons.  It  follows  from 
(9c,  d)  that 

V2  h  »  and  from  the  relation  7^  |v  -  v'l”1  ■  -4ir6^(v  -  v') 

—  a  — 

,  /  m  +  m  \ 

that  V*  h  ■  -4ir  (  - )  f  .  Use  of  these  equations  allows  (9a) 

vo  \  ®a  /  a 


allows  (9a) 


to  be  rewritten  as 
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a  form  useful  for  computational  work  because  the  electron-electron  potential 
h£  does  not  appear. 

I  next  write  the  right  side  of  (10)  in  coordinate  form,  beginning  with 

3 

the  electron-proton  terms.  Under  the  stated  idealizations  fp  »  n  6  ( v) ,  so 
that  (9c)  implies  hp  -  /l  +  £  <v  ,  aiMj  (94)  implies  g  »  nv. 

\  p/ 

Using  (5)  and  (6)  I  now  find  Chat 


— 2.  -  _  JL  e 

3  v  2  v 
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With  f  *  f  (z,  9,  v)  I  also  find  that 
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3ge  32ge 

with  similar  expressions  for  and 


Ic  remains  only  to  perform  the  contractions  indicated  in  (10). 

Neglecting  the  term  4ir  —  f  f  and  changing  variables  from  6  to  u  ■ 

®p  P 

1  find  the  electron  proton-scattering  term  to  be 


For  the  electron-electron  terms  I  obtain 
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(15a) 


(15b) 


(15c) 


(15d) 


(15e) 


p  ■  4  ”  f  ' 


(15f) 


For  notational  convenience  I  have  dropped  Che  subscript  from  gg  in  (15a  -  e) . 

Now  consider  the  inelastic  collision  term  in  (2).  Let  the  imaginary 
trace  element  have  two  bound  levels  with  energy  difference  E^2  and  number 
densities  N^,  and  let  a^2(v)  denote  the  inelastic  cross  section  (assumed 
to  be  isotropic).  For  simplicity  I  take  •  0,  in  which  case  (cf.  Bernstein 
1979,  Shoub  1977) 


6_f 

St 


1 

_2 

rl 

-  N 

in el.  | 

T  °12(V)  I 

I  dp*  f(p',  v)  -  v  o12(v)  f(p,  v) 
-1 

(16) 


where,  by  virtue  of  energy  conservation. 


1-2  1  2 
ymv  «ymv  +  E12 


It  follows  from  (16)  that  the  rate  per  unit  volume  at  which  energy  is  lost 
from  electrons  due  to  inelastic  collisions  is 


-  Jd3  v  (fm  v2)  f£  ^  -  2*1^  E12  J  dv  v3  o^iv)  J 


duf(u,  v)  , 
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-1 


(17) 


(2  M 1/2 

where  v^2  *  y  “  J  •  is  determined  by  equating  the  right  side  of  (17) 
to  the  rate  per  unit  volume  at  which  energy  is  radiated  by  the  plasma. 
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2 

As  mentioned,  1  take  the  latter  rate  to  be  n  $  (T) ,  where  n  is  the  electron 
density  and  ;{>  (T)  is  the  function  given  by  McWhirter  et  al.  (1975).  I  have 
arbitrarily  chosen 


0 ;  v  ^  v 


12 


°12(V) 


*»„  i  v  »  vu 


(18a) 


_q 

where  a  ■  5.29  x  10  cm,  and  we  have  let  E. ,  ■  E,  _  (z)  ■  4kT  (*) ,  where 

O  Li  Li  o 

T  (z)  is  defined  below.  (18b) 

o 

Although  the  above  treatment  of  the  inelastic  collision  terms  is 

schematic,  it  adequately  accounts  for  the  effect  of  radiative  losses 

on  the  electron  temperature  profile.  On  the  other  hand,  it  is 

inadequate  for  treating  possible  distortion  of  the  electron  energy 

spectrum  by  Imbalanced  inelastic  collisions,  and  it  is  therefore 

pertinent  to  ask  if  such  distortion  will  occur  in  the  low  TR.  The 

answer,  most  likely,  is  that  it  will  not.  Any  tendency  of  imbalanced 

inelastic  collisions  to  form  bumps  or  dips  will  be  overpowered  by  the 

smoothing,  Maxwellian-producing  action  of  electron-electron  collisions, 

provided  the  collision  frequency  of  the  latter  exceeds  the  Inelastic 

collision  frequency.  For  thermal  hydrogen  plasma  this  is  the  case 
^  2 

provided  —  <  10  (Shoub  1977),  where  n„  is  the  neutral  hydrogen 
n  n 

e 

density.  In  the  upper  chromosphere,  n^  <  n#  (n^  <  20  n^  in  the  region 

where  hydrogen  Lyman-a  formed),  and  electron-trace  ion  inelastic 

collision  frequency  is  very  small  relative  to  the  elastic  collision 

frequency.  Thus,  except  in  the  event  that  non-thermal  enhancements  of 

3  4 

the  hydrogen  1-2  excitation  rate  (§  IV)  are  as  large  as  factors  of  10  -  10  , 

distortion  of  the  electron  energy  spectrum  by  Inelastic  collisions  is 
unlikely  to  be  important. 
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Dimensionless  Variables 

In  going  from  the  chromosphere  (T  >  10^K)  to  corona  (T  -  2  x  106K) 
the  electron  thermal  velocity  increases  by  a  factor  of  roughly  14.  For 
computational  purposes  it  is  therefore  convenient  to  work  with  dimensionless 
variables.  Toward  this  end  I  introduce  reference  temperature,  density  and 

thermal  velocity  scales  denoted  T  (z),  n  (z)  and  v  .(z)  •  (2kT  /m)1^2 

oo  tn  o 

respectively,  and  choose 

no(*>  Tq(z)  ■  P  •  constant,  (19) 

as  is  appropriate  for  isobaric  slabs.  In  principle  T  (z)  may  be  any  monotonic 

o 

function  satisfying  T  (z  ■  0)  ■  T  and  T  (z  »  L)  »T.  ,  where  T  and  T.  are 

o  co  n  c  n 

temperatures  characterizing  the  incoming  distribution  functions  at  the 
boundaries.  In  practice  it  has  been  found  most  economical  to  have 
Tq(z)  conform  as  closely  as  possible  to  the  actual  temperature  structure  of 
the  atmosphere.  I  therefore  choose  Tfl(z)  as  the  solution  of  the  macroscopic 
electron  energy  equation 


—  (— 

dz  \  in  A 


5/2 

Tq(z) 


dT  \  , 

_ o\ .  2 

— — )+  n 
dz/  c 


*(T  )  - 
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(20) 


T  (z  ■  0)  ■  T 
o  c 

V*  *  L)  ■  th 


Here  K  •  1.87  x  10  5  erg/  (°K^-cm-sec)  is  the  constant  in  the  Spitzer-Harm 
thermal  conductivity,  L  is  the  slab  thickness,  A  ■  3(kT  )2  (4irP)-1^2, 

O  € 

and  $(Tq)  is  McWhirter  et  al.'s  heating  function,  modified  to  allow  for  the 
partial  ionization  of  hydrogen  as  described  in  their  paper.  I  represent 
i  analytically  using  a  piecewlse-linear  fit  (log  0  vs.  log  Tq)  to  their 
graph.  Equation  (20)  is  solved  numerically  using  finite  differences  in 
conjunction  with  generalized  Xewton-Raphson  Iteration. 


I  note  that  the  T  (z)  determined  via(20)  will  differ  from  the  T(z), 
o 

the  temperature  determined  from  the  Fokker -Planck  equation,  only  to  the 
extent  that  the  true  heat  flux  differs  from  that  predicted  by  the  Spitzer-Harm 
formula.  Even  the  latter  T(z)  is  not  self-consistent  however,  due  to 
the  dependence  of  the  radiative  losses,  which  I  am  taking  as  given,  on  the 
tail  of  the  electron  distribution  function.  Except  for  a  few  brief  remarks 
in  SlV,  I  do  not  address  this  latter  problem  here. 

Now  consider  my  choices  for  the  parameters  L,  Tc  and 
T^.  Ideally  one  would  locate  the  upper  boundary  at  the  point  of  maximum 
temperature  in  the  corona  (i.e.,  at  R  ^  3-4  R  ) ,  thereby  accounting  for  the 
entire  reservoir  of  hot  plasma  capable  of  supplying  fast  electrons  to  the 
underlying  TR.  And  the  lower  boundary  would  ideally  be  located 

deep  in  chromosphere  in  order  to  insure  that  all  electrons  thermallze 
within  the  slab.  Unfortunately,  to  model  this  expanse  of  atmosphere 
correctly  requires  consideration  of  several  effects  -  pressure  gradients, 
spherical  geometry,  self-consistent  ionization  equilibrium  calculations  for 
hydrogen  and  the  other  major  electron  donors  -  which  were  judged  best 
omitted  from  this  exploratory  calculation.  Thus  I  have  chosen  the  thickest 
slab  for  which  constant  pressure  remains  a  reasonable  assumption 

4 

(L  *  5  x  10  km) ,  and  a  lower  boundary  temperature  high  enough  that  hydrogen 

remains  almost  fully  ionized  within  the  slab  (Tc  ■  8100K) .  I  set  the  upper 

boundary  temperature  equal  to  the  average  coronal  value  of  2  x  10^K. 

Numerical  results  are  reported  for  these  values  of  T  ,  T.  and  L,  and  for 

c  n 

several  pressures. 

Returning  to  our  discussion  of  the  Fokker-Planck  equation,  I  introduce 
the  dimensionless  variable: 
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(21a) 


(21b) 


(21c) 


(2  Id) 


(21c) 

(21f) 

(21e) 


To  facilitate  comparisons,  I  note  that  a  ■  y  B^,  p  *  -  j  Ag  ■  — ,  where 

c 

Ag  and  B^,  are  the  parameters  appearing  in  Spitzer  and  Harm's  (1953)  paper 

and  E  is  Dreicer's  (1959)  critical  field.  The  parameter  n  enters  the 
c 

PPE  as  a  factor  of  the.  inelastic  collision  terms.  Its  magnitude  can  be 

estimated  by  equating  the  right  side  of  (17)  (evaluated  with  f  a 

2  2 

Maxwellian,  0,.  *  ira  and  E.  „  ■  4  kT  (z))  to  n  $(T  )  and  making  use  of 
1Z  o  12  o  o  o 

-20  -1/2  -1  -1 
the  estimate  (McWhirter  et  al.):  <KTq)  *  5  x  10  Tq  erg  cm  sec  , 

4 

which  is  claimed  to  be  accurate  to  within  a  factor  of  three  for  1.5  x  10 

K  <  T  <  107K.  The  result  is  n  -  1.9  x  10-6. 

—  o  — 

The  electron  thermodynamic  variables  are  related  to  the  reference 
scales  Tq,  nQ  and  through  the  following  expressions  (written  for  zero 
average  electron  velocity) : 


density: 
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/  *  / 
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e2  ♦  (p.  t.  T>dc 


average  velocity: 
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internal  energy: 
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m 


21 


pressure  Censor: 
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heat  flux: 
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Finally,  I  note  that  if  the  distribution  function  is  locally  Maxwellian  with 
density  no(z)  and  temperature  T^(z) ,  its  dimensionless  form  is 


k*Cp.  C,  t)  ■  e”^  • 


i*  satisfies  the  normalization 
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Xhen  written  in  terms  of  the  above  dimensionless  variables  the 
okker-Planck  equation  becomes 
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dp’  <Ku’,  c  +  2,  t)  -  5H(C  -  2)  *<11,  C,  t)' 


The  coefficients  in  (24b)  are  given  in  (15a  -  e)  provided  v  and  g  in  latter 

formulas  are  replaced  by  their  dimensionless  counterparts  £  and  Jl, 

1  3 

respectively;  e.g.,  a  ■  -r- — In  obtaining  (24d)  from  (16)  I  have  used 

95 

(18a,  b)  and  have  introduced  the  notation  H(x)  for  the  unit  step  function. 

A  formula  for  convenient  numerical  evaluation  of  the  potential  £  can 
be  obtained  by  manipulating  the  generating  function  of  the  Legendre  poly¬ 
nomials  (Morse  &  Feshbach  1953,  p.  597)  to  find  that 


(25) 


U  -  C\  <e.  «*>  Pn<cos  Y)  . 

n*0 

where  y  is  the  angle  between  £  and  £.'  and 


Here  £  is  the  larger  and  C<  the  smaller  of  (  and  £'• 


On  substituting  (25)  into  (21d)  and  using  the  spherical  harmonic 

addition  theorem  to  express  P  (cos  y)  in  terms  of  P  (y)  and  P  (y') 

n  n  n 

I  obtain 
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In  practice,  series  (27a)  converges  rapidly  because  £  is  determined  mainly 
by  che  behavior  of  $  at  thermal  velocities  (see  21d) ,  where  <t>  remains  nearly 
isotropic.  An  added  virtue  of  this  formulation  is  that  various  u  and 
£ -derivatives  of  l,  needed  for  evaluating  the  coefficients  (15a  -  f ) ,  can 
be  obtained  analytically  from  the  foregoing  formulas,  thereby  avoiding 
errors  inherent  to  numerical  differentiation. 


Auxiliary  Conditions 

At  the  boundaries  I  choose 

<Ku,  £,  t  -  0)  =  <J>*  (y,  £)  -  4>*(£)  ,  0<u<l 


<t>(y,  £,  t  -  t  )  =  <J>~  (u,  £)  “  $*(£)  ,  -1£U<0 

max  d 


(28) 


where  <t>*  is  a  local  Maxwellian  defined  in  (23a) .  This  choice  can  be  defended 
on  two  grounds  beyond  simplicity.  First,  due  to  the  boundaries  being  located 
in  regions  of  weak  gradients  (top)  and  high  density  (bottom),  I  expect  the 
actual  Incoming  distributions  to  be  close  to  local  Maxwellians.  Second, 

I  find  that  a  few  mean  free  paths  away  from  the  boundaries  the  computed 
solution  is  relatively  insensitive  choice  of  provided  they  do  not 
differ  significantly  from  $*.  For  example,  replacing  <t>*  in  (28)  by  the 
appropriate  Spltzer-Harm  distributions  has  little  effect  on  the  computed 
result  in  the  slab  interior.  The  role  of  boundary  conditions  is  elucidated 
in  paper  II. 

Consider  now  the  task  of  prescribing  constraints  on  the  solution  in 
velocity  space.  To  begin,  observe  that  since  an  electron  at  rest  has  no 
associated  direction  of  flight,  must  be  single  valued  at  £  -  0,  and  hence 

111 

3u  i 


■  0  for  all  ;i,  x. 


(29) 
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Two  additional  constraints  one  wishes  to  impose  on  the  solution  are 
that  it  remain  nonsingular  (particularly  at  £  *  0;  see  Eq.  43)  and  that 
it  decay  at  large  C.  Unfortunately,  there  appears  to  be  no  unique  way 
to  implement  these  constraints  mathematically.  Although  the  forms 
chosen  below  are  physically  motivated,  they  are  not  unique;  others  could 
probably  serve  equally  well.  Fortunately,  numerical  experiments  show  that 
the  solutions  are  insensitive  to  the  particular  forms  chosen  provided  they 
are  sensible  expressions  of  the  above  physical  ideas. 

Proceeding,  I  note  that  (24a)  can  be  recast  as 


■<  {£-!  + 


where  J  is  the  total  velocity-space  flux  density: 


J  -  J  +J  +J,,-6e$ 
—  — ce  — ep  — mel .  z 


It  can  be  verified  that 


4s 


m  .  3  /n2  .  \  1 

2m  3S  V£  Zs)  ~  2 
6  _  w 


32i 

l  3  s  d± 

—  •  -rf  ;  s  •  e,p 


2  3S.3S  ai 


(30a) 


(30b) 


(30c) 
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In  the  following  I  neglect  the  contribution  of  to  J,  because  it  both 

encumbers  the  argument  and  is  physically  unimportant.  Now  let  V  denote  a 


sphere  in  velocity  space  centered  at  the  origin  with  radius  £ 


Integrate 


(30a)  over  V,  use  the  divergence  theorem,  neglect  J 

and  note  that  B(t)  has  presumably  been  chosen  so  that 
(zero  current  condition)  to  obtain 


' )  Y  ) 


d  ££♦(£»  t  )  -  o 


where  T -dependencies  have  been  suppressed.  In  order  to  obtain  a  constraint  on 

the  solution  at  £  I  assume  that  at  each  T  there  are  no  electrons  incident 
max 

on  the  spherical  surface  £  »  ^  from  its  exterior;  i.e.,  I  assume 

—TOT (p  *  ^max^  *  -  0  for  a11  w*  wher*  2T0T  is  the  bracketed  quantity  in  (31). 

It  then  follows  that  (31)  can  be  satisfied  only  if  the  Integrand  vanishes^: 


all  p,t. 


(32) 


The  stipulation  that  no  electrons  be  incident  on  the  surface  £  »  £ 

r  max 

from  remote  regions  of  velocity  space  is  meaningful  provided  the  actual 

number  of  electrons  arriving  at  location  x  with  energy  exceeding 
12  2 

-r  m  v  (t)  *  kT  (x)  £  is  negligibly  small.  A  sufficient  condition  for 
2  max  o  Tnax  66  J 

1  The  term  proportional  to  a  in  Eq.  (32)  was  inadvertently  omitted  in 
calculating  the  results  reported  here.  I  do  not  believe  that  this  is 
a  serious  error. 
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this  to  hold  is  £2  S3CT./T  )**  (say),  for  then  there  are  essentially  no  electrons 
max  h  c 

1  2 

anywhere  in  the  atmosphere  with  energy  exceeding  y  m  vmav  (T  »  0) . 

6  3 

For  1,  ■  2  x  10  K  and  T  ■  (8.1  x  10  K)  this  requires  £  _  >  47.  For  the  present 

n  c  max 

purpose  of  determining  whether  significant  departures  from  LTE  occur,  this  estlmat 
proves  unnecessarily  large,  for  I  find  that  the  computed  solution  for  £  <  £  -  1 

(say)  is  largely  Independent  of  £  provided  £  >5.  However,  the  value 

PMC 

of  £mav.  finally  adop  ted,  namely  £ ■  6,  turns  out  to  be  too  small  to 
accurately  determine  the  electron  hea  t  flux  in  the  lower  TR. 

Again  let  V  denote  a  spherical  volume  in  velocity  space  centered  at 
the  origin,  but  now  with  radius  £  -  e  «  1.  By  Integrating  (30a)  over  V, 
taking  the  limit  e  ■+  0,  and  demanding  that  <p  and  its  derivatives  remain 
bounded  as  £  -*■  0,  I  deduce  that 


lim 

e-*0 


*  [^ee^’  c*  “  6  ♦Cw, 


0  , 


(33) 


which  states  that  the  net  efflux  of  electrons  from  the  origin  is  zero.  To 
obtain  a  constraint  on  the  solution  I  replace  (33)  by  the  stronger  require¬ 
ment  that  electrons  neither  enter  nor  leave  velocity  space  through  the 
origin;  i.e., 


iS  •  { 4e<“-  c>  - s  *<"• 


(34) 


0  ,  all  v»,t. 


.  % 


By  expanding  expression  (30c)  for  J  ,  Eq.  (32)  can  be  brought  into 


the  form 


sft+s£  +  cb*] 


all  u ,  T, 


(35a) 


where 


j2i  . 

7?  ■  H 


* 

i  -  u2  3  at  t 

“  ^2  ay  n  "  C 


(  oE2\  a  1  i  (t2  h\  .11 

*  2w  \s  +  2  /  “  ae  ^2  ac  acj  ^2  ay 


(l  -  y  ) 


(35b, c) 


(35d) 


Equation  (34)  can  be  simplified  by  using  (27a)  to  evaluate  !(  y,  v  ) 
and  its  derivatives  as  5  0.  The  result  is 


+  *b  ♦ 


all  y,  t 


(36a) 


where 


At1  _L  r1  T>  ,  x 
3  *0  15  12  ?2(y) 


(36b) 


{•-!*?} 


(36c) 


Here  I  have  used  Che  nocation 


CO 

-  J  e  -ao  dc. 


where  the  ere  defined  in  (27c) . 

Equations  (28),  (29),  (35a-d)^  end  (36e-d)  ere  the  auxiliary  conditions 
used  in  this  study.  Mote  that  (35a)  and  (36a)  are  satisfied  identically 
by  $  -  0*  if  a  -  8  •  0  and  if  the  coefficients  A^>  ...»  are  evaluated 
for  Maxwellian  field  particles. 


As  a  final  remark  on  auxiliary  conditions,  I  note  that  Kileen  and  Mark 


-1 


(1970)  have  shown  that  the  conditions  (6  -  cos  u) 


H  (•  -  f  .  t  -  0)  -  0 

If  if  o 


are  necessary  for  <J>  to  be  azimuthally  symmetric  about  e^.  They  also  argue 
that  because  |^-  »  -  sin  8  (37b)  does  not  constraint  at  u  *  +1,  and 

consequently  it  is  preferable  to  use  8  rather  than  y  as  the  Independent 
angular  variable.  1  note,  however,  that  (37b)  will  be  satisfied  provided 
remains  bounded  at  y  ■  +1;  this  latter  condition  is  therefore  the 
appropriate  one  when  using  y  as  independent  variable.  I  also  observe  that 
(36a)  reduces  to  (37a)  at  y  ■  0. 


1 


See  footnote  p.  26. 
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Maxwellian  Field  Particle  Approximation 

Ic  will  be  useful  lacer  on  Co  have  explicit  expressions  for  Che  co¬ 
efficients  in  (24b)  when  Che  potential  l  is  evaluated  using  $  ■  <j>*.  This 
has  been  referred  to  earlier  as  the  Maxwellian  field  particle  approximation 
(MFPA).  From  (21d)  and  (23a)  I  find  by  straightforward  integration  that 


>.*«)•£  A3  li-i1  !♦*«’> 


*  (5  +  A) 


erf  (O  +  |- 


Using  this  result  and  equations  (15a-e)  the  coefficients  are  readily 
evaluated.  I  obtain 


erf(g] 

o.3 


* 

_  !_ 
<2 


(39a) 


(39b) 


1-R 


q(C)  ;  q(0 


•(-£) 


erf (^)  + 


(39c) 


d*  - 

r 


(39d) 


e ,v  ■ 


(39e) 


P*  -2  i, 


( 39  f ) 
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Neglecting  terms  of  order  e  %  in  the  above  coefficients  and  equating 
the  proton  and  reference  densities  leads  to  the  following  large-velocity 
form  of  the  Fokker-Planck  equation  in  the  MFPA: 

•<!«-!(»••«)) -£$•(*-;?)» 


To  the  same  approximation  the  auxiliary  condition  (35a)  is 


H + [2  +  **2  <2e  +  a52)]  *} 


For  £  «  1,  the  MFPA  form  of  the  Fokker-Planck  equation  is 


(  „  L  ii+  i  -  "2  ill  .  2  1  +  X  ii  +  2* 

}37-2i5*+t  3tj)  -S("  U  +  —  3  „*  3C  X  2* 


2  /  1  A  S*  1  \  »  /,  2,  Q  6* 

7-  bF  +  T?  (1-“  >  jj  +  sT 

**  7  L  J  inel. 


Lastly,  the  MFPA  form  of  the  auxiliary  condition  (36a)  is 


[— i_  —  +  2w6*  1 
[3  /H  H  J 


Note  that  if 
satisfied  by  0  * 


a  “  6  ■ 


■  0,  nLl  of  the  above  equations  arc 


I  ine l. 
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III.  Numerical  Method 

By  now  it  is  perhaps  apparent  that  the  central  difficulty  in  this  problem 
lies  in  devising  a  means  of  extracting  reliable  information  from  Eq.  (24a). 

The  approach  outlined  here,  a  doubly  iterative,  fully  implicit,  fourth-order 
finite  difference  scheme,  evolved  out  of  my  experience  with  several  simpler  but 
ultimately  unsuccessful  approaches. 

As  mentioned  earlier,  I  treat  the  functional  dependence  of  the  coeffi¬ 
cients  in  the  differential  equation  and  auxiliary  conditions  on  the  distribution 
function  by  iteration,  starting  from  the  approximation  $  ■  4>*  everywhere  (MFPA) . 
This  is  the  first  of  the  two  levels  of  iteration  referred  to  above.  Its 
effectiveness,  is  discussed  at  the  end  of  this  section. 

My  approach  to  the  linear  problem  has  been  motivated  by  the  following 
analogy.  Consider  (24a)  with  the  inelastic  terms  omitted  for  the  moment. 

This  equation  is  of  the  general  form 


“S  If  ■  <A  +  B  "'./'V'1,*0'  W) 

=  L* 

where  the  coefficients  depend  on  y.  5,  and  t.  I  find  empirically  that  A  >  0, 

C  ;  0,  and  4AC  -  B4"  >  0  for  all  y,  5,  x,  except  at  y  *  cl  where  C  *  B  ■  0. 

These  inequalities  can  be  verified  explicitly  in  the  MFPA;  their  relevance  is 
made  clear  below.  Let  ;  =*  :  +  ,  L  «  L+  for  >  0  and  i  *  c  ,  I.  =  L  for  y  <  0 
and  rewrite  (45)  as  the  pair  of  equations 


+ 

♦ 


(46a) 


+ 

o 


(.. 


0)  -  $b  (y.o 


01  u  i  1 


33 


and 


u£  I"-  •  L  ^  (46b) 

♦  v)  *  $v(u,e)  »  -l<u<0 

max  D 

Here  I  have  appended  boundary  conditions  (28)  in  order  to  bring  out  the 
similarity  of  these  equations  to  the  generalized  two-dimensional  diffusion 
equation 


it  *  (A’  3xx  +  B’  axy  +  C'  3yy  +  D’  *x  +  E'  3y  +  «’>*  8  <*>*>  eR‘  t>0 


A*>0  ,  C>0  ,  4A’C'-B2>0  ,  all  (x,y)  eR 


*(x,y,t-0)  -  ¥o(x,y)  , 

where  t  is  time,  R  is  a  region  of  the  x-y  plane,  and  where  now  the  inequalities 
assure  that  the  equation  is  parabolic  (see,  e.g.  Carrier  and  Pearson  1976,  p.  266). 
This  latter  statement  implies  that  Eq.  (47)  can  be  stably  integrated  forward,  but 
not  backward  in  time.  Intuitively  speaking,  this  is  because  the  differential  opera¬ 
tor  on  the  right  in  (47)  is  a  smoothing  operator;  its  action  on  a  spatially  local¬ 
ized  function  produces  a  less  localized  one.  Hence  distinct  initial  configurations, 
vQ(x,y)t  evolve  into  a  similar  ones  as  time  progresses.  Any  attempt  to  reverse 
this  evolution  is  therefore  unstable. 


In  our  case  the  spatial  variable  ds  ■  4^  plays  a  role  analogous  to  time 
in  (47),  while  the  velocity  variables  y  and  5  are  analogous  to  x  and  y. 

Equation  (46a)  is  integrated  forward  in  t,  starting  from  the  boundary  condition 
at  t  «  0;  Eq.  (46b),  on  the  other  hand,  is  integrated  backwards  in  t,  starting 
from  the  boundary  condition  at  Tmav.  An  essential  complication  exists,  however, 
in  that  the  two  equations  are  coupled  through  continuity  requirements  on  $  and 
its  y-derivatives  at  u  »  0.  Fortunately,  this  complication  is  readily  handled 
at  the  difference-equation  level. 

Consider  a  centered,  three-point  difference  approximation 

32<t” 

to  V  near  y  *  0.  Choose  an  equally  spaced  mesh  on  [-1,1]  with  mesh  width 
3y 

Ay  -  2.0/ (NA-1) ,  where  NA  is  an  even  integer.  Let  kQ  -  NA/2  and  label  points 
so  that  “1  ^  y^»  ....  <  0  <  •  •  • »  <  V^jA  £  1.  (Note  that  y  ■  0  is  not 

a  mesh  point.)  Then 


32<t>~ 

3y2 


y-y, 


»k  -l  '  2+k  +  *k  +1 
o _ o  o  A 

(Ay)2 


where  means  :(y^).  Other  y-derivatives  at  mesh  points  adjacent  to  u  =*  0  are 
handled  in  a  similar  way. 

It  is  now  not  difficult  to  see  that  if  the  y  and  ^-derivatives  in  Eqs.  (46a, 
b)  and  side  conditions  (38c,  d)  are  replaced  by  difference  quotients,  with 
u-derivatives  at  mesh  points  near  y  ■  0  treated  as  shown  above  and  those  near 
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u  *  =1  approximated  by  one-sided  differences,  the  resulting  difference  equations 
can  be  cast  in  the  matrix  form 


3£  -  -  + 

—  -  M  £  +  N  £ 


-l<y<0 


(48a) 


3£  +  +  +  - 

-f7  -  M  i  +N  £ 


0<y<l 


(48b) 


where  I  have  divided  (46a, b)  by  the  product  y£  ^  0.  Here  $“(t)  are  vectors 

of  length  kQ  x  NV,  where  NV  is  the  number  of  uniformly  spaced  velocity  points. 

Elements  of  4>~  are  ordered  so  that  all  angle  points  at  a  given  velocity  are 

adjacent;  i.e.,  for  1  £  k  _<  kQ,  the  (k  +  (j  -  1)  k£3)C*1  element  of  £  is 

p(y.  ,  v.,  t),  and  for  k  +  1  <  k  <  NA,  the  (k  -  k  +  (j  -  1)  k  )th  element  of 
kj  o  —  —  o  ■*  o 

£  (t)  is  v j ,  t) .  M~  and  11“  are  square  matrices  of  order  kQ  x  NV,  and,  for 

three-point  velocity  differences,  are  block-tridiagonal  in  form. 

Equations  (48a)  and  (48b)  are  now  solved  by  iteration.  One  starts  with 

_  _  * 

an  initial  estimate  of  £  ,  say  £  *  £  ,  uses  this  estimate  to  evaluate  the  term 
N+£  appearing  on  the  right  in  (48b)  for  all  t,  and  then  integrates  the  result¬ 
ing  inhomogeneous  equation  forward  in  t,  starting  from  t  *  0,  where  i+  *  *  , 

—  — o 

and  ending  at  :  ■  :  The  $+  so  obtained  is  then  used  to  evaluate  the  term 

max.  — 

appearing  on  the  right  in  (48a) ,  and  the  resulting  inhomogeneous  equation 
is  integrated  backward  from  Tmax  to  t*0.  The  process  is  then  iterated  to 
convergence. 

Figure  1  illustrates  how  this  scheme  works  in  practice.  There  I  plot 
successive  approximations  to  $(y)  at  fixed  £(■  5.48)  and  t(«  147;  T  (t)  ■ 

5.28  x  10^K)  for  a  TR  model  with  P  •  P /2.  The  numbers  labeling  curves  denote 
-iteration  number.  Thus,  "0"  is  the  initial  estimate  c*(E  *  5.48),  "1"  is  the 
result  of  the  first  forward  pass,  "2"  is  tie  result  of  the  first  backward  pass. 


and  so  on.  The  curve  labeled  "15"  is  the  result  of  seven-and-a-half  full 
iterations;  as  seen,  it  is  remarkably  smooth.  Behavior  at  other  depths  and 
other  velocities  is  similar.  A  summary  of  my  experience  with  this 
iterative  process  is  the  following.  The  first  few  passes  produce  large  rela¬ 
tive  changes  in  the  solution.  After  ten  or  so  iterations  the  largest  relative 
change  between  iterates  is  generally  a  few  percent;  after  twenty  it  is  a  few 
tenths  of  a  percent.  Thereafter  the  rate  of  convergence  is  slower. 

In  practice  the  iteration  is  stopped  when  the  largest  relative  correction 
to  ':he  previous  iterate  is  a  few  tenths  of  a  percent,  as  this  is  estimated  to 
be  the  levef  of  the  overall  truncation  error. 

Consider,  now  what  turns  out  to  be  the  crux  of  the  numerical  problem, 
namely  choice  of  a  spatial-differencing  scheme.  By  analogy  with  the  situation 
for  time-differencing  Eq.  (47),  there  are  three  classes  of  method  available  (see 
e.g.,  Richtmyer  and  Morton,  1967):  explicit,  alternating  direction  implicit 
(hereafter  ADI) ,  and  fully  implicit  differencing.  Note  that  the  computa¬ 
tional  cost  per  step  of  these  methods  increases  sharply  in  order  listed.  Pertinent 
aspects  of  the  problei  to  be  considered  in  choosing  among  these  are:  a)  our 

interest  in  long  "integration  times"  -  x  can  be  as  large  as  10^;  b)  the 

max 

_  4 

necessity  or  storing  the  solution,  typically  1.5  x  10  numbers  per  depth;  and 

c)  the  strong  velocity  dependence  of  the  coefficients  in  the  differential  equa¬ 
tion;  in  particular,  their  singular  behavior  as  5  0.  Factors  (a)  and  (b) 

imply  that  methods  requiring  a  small  step  size,  although  perhaps  competitive 
from  the  standpoint  of  computation  time,  may  be  unsuitable  due  to  storage 
limitations.  Indeed,  due  to  (c) ,  this  situation  actually  arises  for  both 
explicit  and  ADI  difference  schemes  (see  below). 
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In  practice,  then,  only  fully  implicit  schemes  have  proven  effective, 
essentially  because  only  these  permit  use  of  step  sizes  large  enough  for  the  cal¬ 
culation  to  be  feasible.  The  situation  is  summarized  by  the  following  estimates 
of  the  maximum  local  step  size  At  allowed  by  the  different  methods: 


explicit 

ADI 


(49  a) 

(49b) 


a(t)At  <  1 


fully  implicit. 


(49c) 


Arguments  supporting  (49a, b)  are  given  in  the  Appendix.  Suffice  it  here  to 
say  that  (49a)  arises  from  stability  considerations,  while  (49b)  is  a  necessary 
condition  for  the  difference  between  ADI  and  fully  implicit  approximations  to 
the  spatial  derivative  term  be  uniformly  small.  Since  Ay  %  AS  ^  0.025,  and  in 
light  of  (a)  and  (b)  above,  it  is  clear  that  neither  explicit  nor  ADI  schemes 
are  suited  for  this  problem.  Condition  (49c),  on  the  other  hand, derives  from 
the  underlying  physics  as  can  be  seen  from  the  following  argument. 

As  discussed  in  the  introduction,  electrons  at  location  T  are  collision- 
dominated  for  velocities  such  that  a(r)5  <1.  In  this  regime  the  distorting 

effect  of  gradients  and  electric  field  on  the  distribution  function  is 
mediated  by  velocity-space  derivative  terms  proportional  to  a  and  8  on  the 
left  in  (24a);  the  -r"  term  is  unimportant.  Thus,  at  low  velocities  the  equation 

o  T 

is  local  in  x,  but  is  non-local  in  velocity  space.  On  the  other  hand,  for  large 
velocities  (af;  >>  1)  electrons  traverse  many  temperature  scale  lengths  before 

A 

thermalizing  and  the  equation  becomes  spatially  non-local.  Since  the  —  term 

d  T 

-1/4 

'first  becomes  important  at  ■  a  ,  one  expects  the  largest  allowable  step 
size  to  be  comparable  to  the  mean  free  path  evaluated  at  Ic;  i.e., 


(d  T  I  x 

— — —  J  ,  or  equivalently,  a£i  <_  1.  Since  a  is  quite  small  for 

conditions  of  interest  (Table  1),  it  should  be  possible  to  construct  models 
using  a  modest  number  of  depth  points.  As  mentioned,  only  fully  implicit 
differencing  allows  one  to  realize  this  possibility . 

Unfortunately,  the  computational  cost  per  step  of  fully  implicit  differencing 
is  substantial,  for  at  each  depth  one  must  solve  a  large  system  of  linear 
equations.  For  three-point  difference  approximations  to  ^-derivatives  these 
equations  are  block- tridiagonal  in  form.  Attempts  to  solve  these  equations 
using  iterative  techniques  (successive  block  over-relaxation  and  ADI  iteration) 
proved  unsuccessful,  due  partly,  no  doubt,  to  my  inability  to  find  good  accelera¬ 
tion  parameters  for  either  method,  but  due  also,  in  part,  to  a  more  fundamental 
reason.  Inspection  of  (41)  reveals  that  the  first  y  and  v-derivative  terms  become 
dominant  at  large  velocities.  If  these  are  approximated  using  centered  differ¬ 
ences,  as  is  desirable  from  the  standpoint  of  accuracy,  the  associated  difference 
matrix  becomes  of f-diagonally  dominant,  and  hence  ill-suited  for  treatment  by 
iterative  methods. 

I  was  thus  led  to  adopt  the  direct  method  of  block-triangular-factoriza- 
tion  (block  LU-decomposition) .  In  addition  to  allowing  one  to  obtain  accurate 
soiucions  co  the  difference  equations  in  routine  fashion,  this  approach  has 
the  virtue  of  permitting  many  (up  to  30,  say)  forward-backward  iterations  to 
be  performed  at  a  cost  not  substantially  greater  than  that  of  a  single  pass. 

This  is  because  the  amount  of  computation  required  to  factor  the  original  matrix 
at  a  given  depth  is  far  greater  than  needed  to  solve  the  two  resulting  block- 
triangular  systems.  Disadvantages  of  the  direct  approach,  however,  are  that  it 
requires  extensive  disk  storage  space  and  that,  unless  carefully  coded,  the  cal¬ 
culation  can  become  I/O  bound  during  the  solution  phase  of  execution. 
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Unfortunately,  when  the  above  scheme  was  implemented  using  second-oruer,  cen¬ 
tered  difference  approximations  for  all  derivatives,  and  run  on  the  finest  mesh 
storage  capabilities  would  permit,  it  gave  inaccurate  results.  Subsequent  investi¬ 
gation  showed,  however,  that  accuracy  could  be  improved  substantially  by  going  to 
higher-order  differences.  I  was  thus  led  to  the  following  difference  approxima¬ 
tion  to  (A5)  at  interior  mesh  points: 


u.  C,  a4  Wv  i  -  (a  S4  +  B  S464  +  C  J6  +  D<54  +  E66  +  G 
k  j  r+\Tjk/  \  55  £  u  uv  S  v 


)  M 


uk>0,  i>_5,  kQ<k<NA-3,  3<j<NV-2 


Here  i>  ,  means  <Ku,  ,  £ . ,  t  ),  all  coefficients  are  evaluated  at  the  point 

j  K  K  j  l 

(u,  ,  £  ,  t  ),  and  the  difference  operators  are  defined  as  follows: 

j  l 


£  KJ  -  «4 


ml  ^50  (j1  -  96  9i_1  +  72  $i_2  -  32  ^ ~ 3  +  6  $i-4J 


!i  {7}k|  ■  (U  4S)‘1  [*j-2  -  8  Vl  +  8  Vl  -  ‘j-J 


(51a) 


(51b) 


{‘j'lcl-  (l^c;)2)'1  *  16  -  30  +  16  ;.+1  -  ;3+J]  (51c) 


(60  Au)'1  [-,k. 


3  +  9  *k-2  '  45  Vl  +  45  Vl  *  90  V2  +  Vs 


"“J. 

£  Wk}  ■  [180<i4>2]  [2vs  - 27  \-2  +  270  Vi  -  490  n 

27  V2'  +  2  V3] 


(51d) 


(51e) 


Here  unneeded  subscripts  and  superscripts  have  been  omitted.  Superscripts  on  tne 
r,s  denote  the  order  of  the  associated  truncation  error. 
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Calling  the  right  side  of  (50)  ?  ,  I  write  the  difference  equations 

J  * 

for  grid  points  near  the  lower  boundary  as: 


i  -  2,  3,  yfc  >  0 


(52a) 


i  -  4,  >  0 


(u  ♦*  - 18  »1'1  +  9  *l~2  - 2  ^'3>  ■  »}k 


(52b) 


For  the  y-derivatives  near  end  points  I  use  centered  fourth-order  formulas 

at  y„.  and  one-sided  fourth-order  formulas  at  y...  ,  and  uvt..  ^-derivatives 
NA—  l  NA—  1  NA 

at  C  *  AC  and  C  “AC  are  approximated  by  centered  second-order  formulas,  and 
nicix 

at  C  *  0,  C  by  one-sided  second-order  formulas.  Use  of  these  second-order 
max 

formulas  does  not  seem  to  have  degraded  the  overall  fourth-order  accuracy. 

The  foregoing  difference  equations,  together  with  their  counterparts  for 
_  <  0  and  difference  approximations  to  the  subsidiary  conditions  (42)  and  (44), 
can  be  written  in  matrix  notation  as: 


+  +  •  „+  - 
■lih  m  -  Si  ♦i  . 


y>0,  i«2,  . . . ,  ND 


(53a) 


-  £  -  fi  £  • 


y<0,  i«*ND-l,  ...,  1 


(53b) 


Here  ND  is  the  number  of  depths  and  the  vectors  represent  that  part  of  the 
spatial  derivative  term  depending  on  £7  for  *  <  i  (a  >  0)  or  i  >  i  (-  <  0);  at 
any  stage  of  the  calculation  the  are  therefore  known  vectors.  The  origin  of 
matrices  N*  and  the  ordering  of  elements  of  has  been  discussed  in  the  lines 
following  (48b).  Matrices  are  of  order  kQ  x  NV  (recall  kQ  -  NA/2)  and  are 
block-pentadiagonal  in  form; 

‘Si  ii  fil 

§2  ^2  §2  0 

-3  §3  -3  *3  -3 


'ENV-2 


V  &V  SnV 


The  submatrices  A ,  . . ,  are  square  and  of  order  kQ;  in  general,  they  are  all 
bandmatrices  with  bandwidth  seven.  In  the  MFPA,  however,  off-diagonal  sub- 
matricet  are  diagonal  due  to  the  absence  of  the  mixed  derivative  term  from  (45) 

±  4* 

From  the  defining  relation  -  L"y“  it  is  easily  found-  that 


;  y"i 


It] 
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where  the  submatrices  are  again  square  and  of  order  kQ  and  I  is  Che  identity 
matrix.  Ke  note  that  although  the  original  A ^ are  sparse,  their 

^  A  A 

counterparts  A  ^ (except  Ej  *  E^)  are,  in  general,  dense. 

Several  comments  concerning  my  choice  of  difference  formulas  are  in 
order.  First,  I  note  that  the  submatrices  in  U^and  L^fill  regardless  of  the 
bandwidth  of  the  submatrices  in  M^.  Computation  and  storage  requirements  are 
therefore  largely  independent  of  the  order  of  p-difference  formulas.  In 
addition  to  small  truncation  errors,  the  chosen  seven-point  formulas  afford 
enhanced  coupling  between  positive  and  negative  p-components  of  the  solution. 
Similarly,  the  use  of  five-level  spatial  differences  entails  less  storage  and  com¬ 


putation  than  do  centered,  second-order  formulae,  yet  offer  smaller  truncation 
errors.  Unfortunately,  going  from  three  to  five-point  velocity  differences  costs 
dearly;  for  a  fixed  number  of  velocity  mesh  points,  disk  storage  doubles  and  com¬ 
putation  increases  by  roughly  a  factor  of  2.5.  Nevertheless,  I  find  that  the  increased 
accuracy  of  the  fourth-order  method  more  than  compensate  for  the  added  storage 
requirements  per  point. 

The  author  has  recently  learned  that  it  is  possible  to  difference  (4o) 
in  such  a  way  as  to  obtain  fourth-order  accurate  approximations  to  the  p  and 
'-derivatives  and  yet  have  the  associated  difference  matrices,  M~,  oe  block- 
tridiagonal.  This  procedure  is  being  tested  and  if  proven  effective  will  be 
used  in  future  versions  of  the  present  code.  Also,  newly  developed  iterative 
methods  for  solution  of  non-symmetr ic  sparse  matrices  (Manteuffel  1977,  Kershaw 
1977,  Van  Der  Vorst  1981)  appear  promising. 


r— 


An  overview  of  the  linear  problem  can  be  had  by  viewing  its  global  matrix 


structure: 
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Here  I  is  the  k  x  NV  identity  matrix,  are  vectors  containing  grid  point 
-  o  -  o 

values  of  the  boundary  conditions,  x' s  denote  diagonal  matrices  originating 
from  one-sided  difference  approximations  to  the  spatial  derivative  term,  and 


the  y's  are  band  matrices  originating  from  our  use  of  centered  space  differences 
at  1*2,3  (u  :•  0)  and  at  i  -  XD-1,  ND-2  (u  *■'  0). 


The  forward-backward  iteration  described  earlier  is  now  seen  to  be  the 
obvious  approach  to  the  iterative  solution  of  the  above  matrix  equation. 

The  above  procedure  has  been  programed  ana  tested  on  the  Cray-1  computer 
at  the  National  Center  for  Atmopsheric  Research.  Due  to  the  predominance  of 
matrix  manipulations  the  algorithm  is  Ideally  suited  to  the  vector  processing 
capabilities  of  the  Cray-1.  In  fact,  essentially  all  of  the  "hard"  computa¬ 
tion  is  performed  using  assembly-language  routines  designed  to  utilize  the 
vector  character  of  the  machine  to  the  fullest  possible  extent.  Results 
presented  in  the  following  section  were  obtained  using  71  depths,  60  angles 
and  256  velocities.  For  this  number  of  points  approximately  six  minutes  of 

central  processor  time  is  required  to  setup  and  solve  the  matrix  equation  (54), 

+ 

using  30  forward-backward  passes.  Storing  the  triangular  factors  of  the  Jl“ 

matrices  requires  2  x  (ND-1)  *  NV  *  4  *  k  *  (k  +  7)  %  1.6  x  10  ^  words  of  disk 

o  o 

space.  High  disk-to-central  memory  data  transfer  rates  are  obtained  by  using 
four  independent  data  channels  and  a  double  buffering  scheme. 

The  code  was  tested  in  three  ways:  first,  it  was  verified  that  starting 
from  a  distribution  markedly  different  than  Maxwellian,  the  calculation  would 
converge  to  an  absolute  Maxwellian  in  an  isothermal  and  isobaric  atmosphere 
(.(-)  =  2(t)  *  0)  with  Maxwellian  boundary  conditions.  Second,  it  was  verified 
that  the  Spitzer-Harm  solution  could  be  reproduced  in  a  high  pressure  atmosphere. 
Four  coefficient  iterations  were  required  to  obtain  an  accuracy  of  five  percent 
or  better  for  £  £  3.2.  Finally,  the  accuracy  of  calculation  under  realistic 
conditions  was  checked  by  observing  the  effect  of  judiciously  refining  the 
mesh.  For  this  purpose  a  word-packing  routine  was  used  to  effectively  double 
the  amount  of  disk  storage  available.  Unfortunately,  use  of  this  routine 
-proved  far  too  costly  for  it  to  be  used  for  other  than  testing  purposes. 
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As  a  result  of  Chese  tests  I  judge  the  results  presented  below  to  have  a 

relative  accuracy  of  ten  percent  or  better.  J 

I  turn  now  to  the  promised  discussion  of  the  effectiveness  of  treating 

the  quasi-linear  collision  term  by  iteration  under  realistic  conditions.  Table  2a 

is  a  tabulation,  for  three  coefficient  iterations,  of  the  zeroth  and  first 

angular  moments  of  the  distribution  (cf.  27c),  normalized  to  $*,  at  a  location 

near  the  bottom  of  our  standard  model  (x  ^  27,  Tq(t)  *  19900K,  a(t)  ■  3.8  x 
_2 

10  ).  The  calculation  was  started  from  LTE  so  that  the  first  iterate  corre¬ 

sponds  to  the  MFPA.  The  points  to  note  regarding  Table  2a  are:  a)  the  largest 
relative  change  in  <p°  between  the  first  and  third  Iterates  is  less  than  four  per¬ 
cent  for  all  :£  <  5.5.  b)  the  same  is  true  for  for  £  <  C  <  5.5.  (£  * 

—  c  —  c 

a  ^2.3).  and  c)  fluxuates  widely  between  Iterates  at  thermal  velocities. 
On  the  basis  of  these  and  more  detailed  comparisons,  I  infer  that  <f>°  at  all  £, 
and  the  entire  distribution  function  at  suprathermal  velocities,  can 
be  obtained  with  good  accuracy  in  the  MFPA.  What  cannot  be  obtained  correctly 
in  this  approximation  are  quantities  which  depend  sensitively  on  the  pitch 

angle  distribution  of  thermal  electrons;  for  example,  the  heat  flux.  This  is 

00 

0  r  5  1 

illustrated  by  Table  (2b)  giving  the  dimensionless  heat  flux,  J  5  $  (£,  t)  d£, 

o 

for  each  iteration  together  with  its  Spitzer-Harm  value. 

I  close  this  section  by  mentioning  two  simplifications  which  have  been 
useful.  First,  I  evaluate  the  troublesome  inelastic  term  (24d)  iterativelv. 
during  the  course  of  the  forward-backward  iteration.  Second,  I  determine  the 
polarization  electric  field  required  for  zero  current  flow  using  the  Spitzer- 
Harm  (1953)  result  3*  0.35  a.  This  value  of  3  gives  zero  mean  electron 
velocities  to  well  within  the  accuracy  of  the  calculation. 


IV. 


Numerical  Results 


Results  are  presented  for  two  model  transition  regions,  one  with 

14  -3 

electron  pressure  <n  T  >  ■  6  x  10  cm  K  =  P  and  the  second  at  P  / 2- 

oo  o  o 

1  note  that  although  PQ  is  the  most  frequently  quoted  electron  pressure 
for  the  spatially  averaged  quiet  sun  TR  (e.g.,  Dupree  1972),  different  EUV 
observations  yield  pressures  which  can  differ  from  PQ  by  at  least  a  factor 
of  two  (Witte roe  1978) .  Both  models  span  the  temperature  range  8100  - 

2  x  10^  K  and  have  a  thickness  of  5  x  10^  km.  The  distribution  function 
has  been  calculated  in  the  MFPA  for  the  boundary  conditions  (28)  using 

■  6  and  a  grid  of  256  velocities.  71  depths  and  60  angles, 
a)  Reference  Atmosphere 

A  portion  of  the  reference  atmosphere  obtained  by  solving  Eq.  (20) 

with  <n  T  >  «  P  is  given  in  Table  3,  where  \  is  the  thermal  mean  free 
o  o  o 

path,  qgh  the  Spitzer-HArm  heat  flux,  and  Qrad  the  radiative  losses  per  unit 

volume  per  unit  time  calculated  from  McWhirter  et.  al*s  (1975)  loss  function.  Two 

points  concerning  this  table  should  be  noted.  First,  from  the  near  constancy  of 

qgh  above  roughly  20,000  K  one  can  infer  that  radiative  losses  have  a  negligible 

effect  on  the  temperature  profile  above  this  level.  Consequently  T  (z)  fcr  the 

o 

P  / 2  case  is  very  similar  to  that  for  P  and  values  of  the  parameter  o  are 

o  o 

roughly  double  those  for  the  PQ  case.  Second,  relative  to  empirical  models 
(e.g.,  Vernazza  et  al.  1973,  Gouttebroze  and  Leibacher  1980), 
the  calculated  temperature  gradient  is  probably  too  steep  below  60,000  K, 
and  is  certainly  so  below  30,000  K.  This  discrepancy  illustrates  the  known 
result  (Withbroe  1977,  Pneuman  and  Kopp  1978,  Athav  1981)  that  the  empirical 
temperature  profile  of  the  low  TR  in  either  the  quiet  sun  or  in  active  region 
loops  cannot  be  understood  on  the  basis  of  models  whose  energy  budget  includes 
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only  ^classical)  thermal  conduction  and  radiation.  However,  in  the  present 

context  this  discrepancy  is  irrelevant  because,  as  mentioned  earlier,  Tq(z) 

serves  only  to  define  the  relationship  between  dimensional  and  dimensionless 

variables;  the  true  temperature  profile  is  determined  self-consistently 

from  the  Fokker-Planck  equation  and,  in  principle,  is  independent  of  Tq(z) . 

This  is  not  wholly  the  case  here  because  radiative  losses  are  (incorrectly) 

evaluated  on  the  T  (z)  scale  and,  moreover ,  are  assumed  to  have  the  same 
o 

functional  dependence  on  temperature  as  in  (electron)  LTE. 


b)  Thermodynamic  Quantities 

Typical  values  of  the  reduced  density,  temperature,  mean  velocity, 

and  the  ratio  of  parallel  to  perpendicular  pressures  are  listed  in  Table  4 

(see  Eqs.  22  a-f) .  Points  of  Interest  are  that:  a)  the  Spitzer-H&rm  value  for 

the  polarization  electric  field  gives  adequately  small  mean  electron  velocities 

in  the  MFPA.  b)  Below  it  will  be  seen  that  the  tail  of  the  distribution  is 

highly  anisotropic  (Fig.  5).  Nevertheless,  the  ratio  p,  (  /p^  remains  within 

a  : ev  percent  of  unity.  Whether  this  remains  true  when” the  calculations  are 

carried  to  larger  velocities  remains  to  be  determined,  c)  The  ratio  T/T  is 

o 

seen  to  be  slightly  larger  than  unity  in  the  lower  atmosphere.  No  quantitative 
significance  should  be  attached  to  this  temperature  rise  however,  because 
the  heat  flux  in  the  present  calculation  is  inaccurate  due  to  my  use  of  the 
MFPA  and  to  truncation  of  the  calculation  at  -  6.  The  latter  shortcoming 

appears  especially  serious  in  light  of  recent  calculations  using  the  Bhatnaghar, 
Gross  and  Krook  equation  (Shoub  1982;  Paper  II),  which  indicate  that  a  substan¬ 
tial  part  of  the  energy  flux  into  the  low  TR  is  carried  by  superthermal 
electrons  with  velocities  10-30  times  local  thermal  velocities.  Thus,  the 
correct  Ckine tic- theoretical)  temperature  profile  of  the  TR  remains  to  be 


dot err ined . 
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Finally,  I  note  that  although  the  calculated  T(z)  is  inaccurate  low 
in  the  atmosphere,  this  does  not  imply  that  the  calculated  distribution 
function  (see  below)  is  grossly  in  error  in  this  region.  In  all  likelihood 
it  is  not,  especially  at  supra  thermal  velocities,  because  the  tail  of  the 
distribution  in  the  low  TR  is  largely  independent  of  the  local  temperature 
structure. 

c)  Distribution  Function 

The  Isotropic  part  of  the  dimensionless  distribution  function,  <f  °  (£,  T  )  * 

° 

ZJ  dp  $(u»  5»  T),  is  shown  at  several  heights  for  the  case  P  ■  Pq  in  Fig.  2, 
and  in  more  detail  in  Table  5.  The  effect  on  d> 0  of  reducing  the  pressure  a 

factor  of  2  is  also  illustrated  in  Table  5.  The  dimensional  distribution  f°  (v,T  ) 

v  ’  o 

3  c 

(no/2ir  v^)  $  (C,  Tq)  is  shown  in  Fig.  3  for  the  case  P  -  p0/2'  ^or.  comParison  a 

local-Maxwellian  distribution  is  shown  as  a  dashed  line  in  Figures  2  and  3. 

The  essential  results  here  are  that,  relative  to  a  local-Maxwellian, 

the  tail  of  the  distribution  is  highly  overpopulated  and  exhibits  only  a  weak 

spatial  (temperature)  dependence  throughout  the  lower  TR  and  chromosphere. 

From  Table  5  it  can  be  seen  that,  with  respect  to  normalization,  the  enhanced 

* 

tail  is  compensated  for  (approximately)  by  a  suppression  of  b  below  $  for 
1.5  <  Z.  <  2.5  For  temperatures  less  than  approximately  10°K,  begins  to 

•k 

exceed  c  at  about  5  *  2.5.  Higher  in  the  atmosphere  <p°  remains  close  to 
* 

4>  out  to  progressively  larger  C.  An  analytical  characterization  of 
<*>°  (',  T^)  is  obtained  in  paper  II. 

Fig.  4  shows  the  1-d‘menaional  distribution 

*z  (;z,T0>*  J1  dM  /  ?2  6  (5z  ■  ♦  (u»  5.  Tc)d$ 

E  *  (E  /E,  E,T  )dE 
z  o 

at  several  heights,  including  the  boundaries.  Here  '  «„•  A 

Maxwellian  is  again  shown  as  a  dashed  line  for  comparison. 


N'oce  che  asymmetry  between  directions  parallel  (C 2>0)  and  anti-parallel 

to  7T.  Unfortunately,  the  structure  in  for  £  >0,  evident  in  Fig.  4  at 

z  z 

Tq  -  19,900  K  and  1.05  x  103  K,  is  artificial — an  artifact  of  my  setting 
£)  ■  0  for  £  >  6  in  the  above  integral  for  <p  .  This  procedure  is 
inaccurate  for  £z  >  3  because  the  pitch  angle  distribution  (Fig.  5,  y  >  0) 
is  such  that  for  £  >  3,  4>(™»  £)  decreases  slowly  with  increasing  £. 

Finally,  I  note  that  <Pz  differs  from  4>*  at  thermal  velocities,  but 
the  differences  are  too  small  to  be  noticeable  in  Fig.  4. 


Fig.  5  shows  the  pitch  angle  dependence  of  the  distribution  for  a 

sequence  of  velocities  at  a  location  near  the  bottom  of  the  TR  for  the  case 

<n  T  >  ■  P  /2.  The  main  features  of  these  plots  are  the  following.  At  low 
o  o  o 

velocities  deviations  from  isotropy  are  small  and,  moreover,  $  varies  linearly 

with  p,  as  predicted  by  classical  theory.  The  flux  of  low  velocity  electrons 

is  directed  upward,  toward  higher  temperatures.  This  upward  drift  has  two 

causes;  natural  (concentration)  diffusion  due  to  a  negative  gradient  in  the 

dlnF  9  c  din  T 

number  of  low  velocity  electrons  (  dg  =  ~dz~  =  (£“■§)  — dz  — ) »  and  forced 

diffusion  driven  by  the  polarization  electric  field.  Evidently,  for  some  £ 

between  1.75  and  2,^~  becomes  positive  and  sufficiently  large  that,  in 

response,  suprathermal  electrons  diffuse  downward  against  the  action  of  the 

field.  This  effect  becomes  more  pronounced  with  increasing  £,  leading  to  the 

strong  anisotropy  evident  in  the  bottom  two  panels  of  Fig.  5.  Note  that  for 

£  >_  3,  i(  y,£  )  decreases  relatively  slowly  with  increasing  y  for  -1 

where  u  (£*3)~  .5  and  y  (£  ■  4) t  .4,  and  then  decreases  more  rapidly  for 
c  o 

y  (£)  <  u  <  1.  That  the  changeover  occurs  at  y  (£)  >  0,  and  not  aty  *  0, 
o  —  o 

is  due  to  pitch  angle  (back)  scattering. 
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Fig.  6  is  a  contour  plot  of  the  distribution  function  at  a  location 
near  the  bottom  of  the  TR  for  the  case  P  =  P  .  From  this  plot  it  is 
clear  tnat  the  propagation  of  electrons  down  the  gradient  is  diffusive 
rather  than  beam- like, in  character. 


d)  Effect  on  Collision  Rates 

From  the  standpoint  of  spectroscopic  diagnostics,  the  main 
consequence  of  having  an  overpopulated  tail  on  the  electron  energy 
spectrum  is  that  it  causes  an  increase  in  electron-ion  collision  rates 
and  an  attendant  alteration  of  an  element's  ionization  balance  and, 
possibly,  its  line  emission.  I  illustrate  these  effects  here  for  magnesium 
and  for  helium,  but  stress  that  the  results  are  preliminary  due  to 
shortcomings  in  both  my  calculation  of  the  distribution  function  and 

ionization  equilibrium  •  (see  below), 
i)  Magnesium 

The  expression  for  the  collisional  ionization  rate  per  ion  is 

C1k  (T)  "  4it  /  V  °1«(V)  f°(v)  dV 
V1k 

00 

•  2  Vth  J  53  °i«  (vth°  4  <«■ T)  « 

where  a,  is  the  ionization  cross  section,  v  is  the  electron-threshold  velocity, 
Ik  •lk 

and  E  ■  v  /v  ..  The  corresponding  rate  based  on  a  Maxwellian,  denoted 
Ik  Ik  th 

*  ,  , 

C^,  is  easily  shown  to  be 

C*  (T)  -  n  C  \/t~  exp  (-E.  /kT)  T(T) 

Ik  o  ak 
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whore  is  a  constant.  ■  1/2  m  and 

r(T)  -  /  (x  +  (E1k  +  xkT)  e'*  dx 

«  ^ 
is  slowly  varying.  In  Table  6  I  list  and  the  ratio  C^/C^  for  the 

first  six  ionization  stages  of  magnesium.  C1k  was  calculated  numerically 

using  the  FP  results  for  and  the  cross  section  from  Allen  (1976  p.  41). 

There  are  no  entries  for  >  5.4  because  4>°  was  calculated  only  out 

* 

to  C  *6.  The  salient  aspects  of  these  results  are:  a)  C,  / C,  *  1 
max  Ik  lie 

4 

when  is  below  a  critical  value  increasing  from  roughly  2.0  at  T  -  10  K 
to  2.5  at  T  «  10^  K.  b)  Once  is  larger  than  this  critical  value 

jig 

C^/C^  increases  exponentially  with  decreasing  temperature,  and  c)  once 
* 

C^/C^  >>  1  at  a  given  temperature,  the  ratio  is  then  strongly  density 
dependent.  These  results  simply  reflect  the  behavior  of  the 

calculated  distribution  function:  (a)  holds  because 

-°( ^ > T0 )  =  i*(£)  for  £  <  Cc  and  Cc  increases  with  T0;  (b)  holds  because  at 
low  temperatures  the  tail  of  the  dimensional,  angle-averaged  distribution, 
and  hence  ,  has  only  a  weak  spatial  dependence  (Fig.  3)  whereas  decrease 
exponentially  with  decreasing  temperature;  (c)  simply  reflects  the  density 
(or  pressure)  dependence  of  the  tail,  as  shown  in  Table  5. 

Figure  7  shows  how  the  enhanced  ionization  rates  affect  the  ionization 
balance  of  Mg.  Processes  included  in  the  ionization  equlibirium  calculation 
were  collisional  ionization  and  radiative  and  dielectronic  recombination. 

The  recombination  coefficients  were  taken  from  Aldrovandi  and  Pequignot  (1973) . 


The  essential  features  of  Fig.  7  are  that:  a)  ionization  stages  normally  found  in 
the  upper  and  middle  TR  remain  unaffected  (see  below) .  b)  The  temperature  range  ov< 
which  lower  ionization  stages  exist  is  generally  broadened  (on  the  low 
temperature  side)  and,  in  some  instances,  shifted  to  lower  temperatures,  and 
c)  a  density  dependence  is  introduced  into  the  ionization  equilibrium. 

Clearly,  electron  non-LTE  effects  significantly  alter  magnesium's 
ionization  equilibrium.  Qualitatively  similar  results  have  been 
found  for  C,  N,  0,  Si,  S,  but  it  is  felt  premature  to  present  extensive 
ionization  equilibrium  data  until  deficiencies  in  the  present 
calculation  are  remedied.  Four  deficiencies  in  particular  may  alter 
the  results  shown  in  Fig.  7.  First,  the  fact  that  the  calculated  tempera¬ 
ture  profile  is  too  steep  at  low  temperatures  implies  that  <J>° departs  from 
o  at  lower  values  of  £  than  is  perhaps  realistic.  This  leads  to 
an  overestimate  of  ionization  rates  for  low  ionization  states.  Second,  ionization 
rates  of  higher  ionization  stages  have  been  underestimated  at  low  temperatures 
because  in  calculating  the  ionization  rates  I  set 
ionization  any  rate  for  which  £  exceeded  5.4  equal  to  zero. 

This  turns  out  to  be  a  poor  approximation  because  the  tail  extends  to  many 
10' s  of  thermal  velocities  in  the  low  TR.  Third,  I  have  not  attempted  to 

account  for  effects  of  the  nonthermal  electron  energy  spectrum  on  dielectronic 
recombination  rates.  However,  because  this  process  involves  inelastic  electron- 
ion  collisions  and,  moreover,  favors  collisions  with  large  excitation  energy 
(e.g.,  Mihalas  1978),  I  suspect  its  rate  will  be  significantly  enhanced  at 
temperatures  below,  say,  3  x  10^  K.  If  it  turns  out  that  the  rate  of  dielectronic 
recombination  exceeds  that  of  radiative  recombination  in  the  low  TR,  the  ionization 
balance  may  be  significantly  different  than  shown  in  Fig.  7.  Fourth,  I  note  that 
provided  the  effect  just  mentioned  does  not  qualitatively  change  present  results,  it 
is  likely  chat  at  low  temperatures  the  situation  will  be  one  in  which  neutral 


hydrogen  and  helium  coexist  with  multiply-charged  trace  elements.  In  this  case 


charge  transfer  recombination  of  trace  ions  in  collisions  with  atomic  hydrogen 
and  helium  (Butler  and  Dalgarno  1980)  may  play  an  important  role  in  the 
ionization  balance. 

From  the  above  discussion  it  is  apparent  that  a  considerable  amount  of 
work  remains  to  be  done  before  reliable  ionization  equilibrium  results  are 
available. 

Now  consider  the  question  of  whether  a  non-thermal  electron  energy 
spectrum  will  cause  enhanced  emission  in  a  given  spectral  line.  Several 
cases  can  be  distinguished.  (1)  If  the  fractional  abundance  curve  of  the 
relevant  ion  has  not  been  broadened  or  shifted  to  lower  temperature 
(e.g..  Mg  VI  in  Fig.  7),  then  enhanced  emission  is  not 

expected.  This  is  because  an  unaltered  abundance  curve  implies  C,  •  C? 

r  lie  lie 

over  the  temperature  range  in  which  the  ion  is  abundant  and,  since 

and  departures  from  a  Maxwellian  spectrum  increase  mono tonically  with  energy, 

* 

this  implies  C12  =  over  the  same  temperature  range,  (ii)  If  the  ion’s 

abundance  curve  has  been  shifted  (e.g..  Mg  IV),  and  if  E^  is  comparable 

to  E^,  then  emission  will  be  enhanced  because  a)  emission  occurs  over 

a  broader,  more  dense  region  than  in  LTE  and  b)  because  Is  likely  to 
★ 

exceed  over  this  range.  Thus,  resonance  lines  of  hydrogen  and 
helium  (see  below)  are  likely  candidates  to  show  enhanced  emission. 

(Lii)  Finally,  if  the  relevant  abundance  curve  has  been  altered  but  E^2 
is  small  compared  to  E^,,  then  although  enhanced  emission  will  probably 
still  occur  due  to  density  and  range  effects,  the  relative  enhancement 
will  be  less  than  in  case  (ii)  because  C^2  will  be  close  to  C^2  over  the 
region  of  formation. 

The  above  remarks  are  clearly  qualitative;  detailed  calculations 
are  required  if  one  is  to  draw  reliable  conclusions. 


54 


r 


s 

s 


11)  Helium 

The  quiet  sun  XUV  spectra  of  helium  has  evoked  much  interest  in 
recent  years  (see  Mango  et  al  1978  for  a  review),  mainly  because  it  is 
poorly  understood.  The  central  issue  has  been  to  identify  the  mechanisms 
responsible  for  ionizing  Hel  and  Hell  at  low  temperatures  (Athay  1975,  p.  298)  and 
for  enhancing  resonance  line  intensities  considerably  in  excess  of  that 
predicted  by  standard  colllsional  excitation.  Although  it  is  now  bellved 
that  a  photoionization-radiative  recombination  mechanism  is  responsible  for 
formation  of  the  resonance  continue  of  both  Hel  and  Hell  (Mango  et  al) , 
the  exciting  mechanism  for  the  Hel  X584  and  Hell  X304  lines  remains  ambiguous. 

In  1975  Jordan  suggested  that  what  was  needed  was  a  method  of  mixing  under¬ 
ionized  helium  with  hotter  electrons,  thereby  enhancing  the  line  emission. 
Subsequently,  Shine  et  al  (1975)  reported  calculations  which  indicated  diffusion 
of  helium  up  the  temperature  gradient  was  one  means  of  achieving  such 
mixing.  I  show  here  that  diffusion  of  fast  electrons  down  the  temperature 
gradient  is  a  second. 

,  * 

In  Table  7  I  list  ionization  rate  ratios,  C..  /C.  ,  for  Hel  and  Hell 

Ik  Ik 

and  excitation  rate  ratios  for  the  Hel  X584  and  Hell  X3D4  transitions.  Both 
actual  and  LTE  rates  were  calculated  numerically  using  cross  sections  from 
Mlhalas  and  Stone  (1968) <  Note  that  the  Hell  rates  are  considerably  enhanced 
even  at  45,800  K,  well  above  the  temperature  at  which  we  expect  thermalizatlon 
to  occur  in  a  more  realistic  chromosphere  model. 

Proper  assessment  of  how  the  collision  rate  enhancements  indicated 
above  affect  helium's  ionization  balance  and  line  emission  requires  radiative 
transfer  calculations  beyond  the  scope  of  this  paper  (see:  Avrett  et  al  1976; 
Milkey  et  al  1973) .  The  following  qualitative  observations  may  be  noted 
however.  First,  it  appears  that  helium's  ionization  balance  will  in  fact  be 
altered  by  the  effects  in  question,  for  the  enhanced  collisional  ionization 
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rates  are  comparable  to  the  nominal  photoionization  rates  for  Hel  and  Hell 
quoted  by  Avrett  et  al.  throughout  the  temperature  range  of  interest. 

Second,  the  observed  brightening  of  the  helium  resonance  lines  (and  other 
EUV  emission  lines)  in  going  from  cell  interior  to  network  region  can  be 
understood,,  at  least  partially,  in  terms  of  the  sun's  large-scale  magnetic 
field  structure;  i.e.,  the  predominantly  horizontal  field  in  cell  Interiors 
impedes  vertical  diffusion  of  electrons.  Third,  the  observed  decrease  in 
brightness  of  the  helium  resonance  lines  in  coronal  holes— which  is  more 
pronounced  than  for  other  ions  (Jordan  1975,  Bohlin  1977) — and  their  relative 
brightening  in  active  regions  can  be  understood  in  terms  of  1)  their  relatively 
large  ratio  of  excitation  to  ionization  energy,  and  li)  differences  in 
physical  conditions  in  these  regions,  i)  implies  the  lines 
have  large  values  of  E^/kT  in  their  region  of  formation  and  hence  are 
particularly  sensitive  to  a  non-thermal  electron  energy  spectrum. 

Regarding  ii)  1  note  that  a  useful  parameter  for  assessing  the  like¬ 
lihood  of  non-Maxwellian  effects  is  collislonal  depth  xL  separating  high  and  low 
temperature  regions.  For  an  isobaric  slab  of  thickness  L,  electron  pressure  F, 

maximum  and  minimum  temperatures  and  Tc,  and  possessing  a  conduction 

7/2 

temperature  profile  (i.e.,  T  linear  in  height),  I  find  that 

L  . 

r  -  f  iz'  s  7L  »  7TTe4  P  InA  ,<A  lnT>\  , 

L  /  T(z*)  A(Th)  (kTh) 3  \  L  / 

o  T  -  T 

where  I  have  assumed  T.  >>  T  so  that  <A  In  T>  z  - = — —  -  1.  Coronal  holes 

h  c  Th 

are  thought  to  be  cool,  low-density  regions  with  relatively  weak  temperature 
gradients,  whereas  active  regions  are  hot,  dense  regions  with  steep  gradients. 
Thus,  T  most  likely  Increases,  and  non-thermal  effects  therefore  decrease, 

Lt 

in  going  from  active  regions  to  holes.  Fourth,  and  perhaps  most  interestingly, 

I  note  that  in  addition  to  its  potential  ability  to  account  for  the  observed 
emission  in  the  '584  and  A304  resonance  lines,  the  present  theory  helps 
explain  a  puzzling  aspect  of  the  helium  spectrum,  namely  the  large  observed 


values  of  the  Lyman  series  intensity  ratios  for  Hell; 

I(n  +  1;  n  >  3)/I(2  ■+  1) .  As  no  ted  by  Zirin  (1975),  it  is  difficult  to 

explain  the  observed  intensity  ratios  via  collislonal  excitation  by 

electrons  with  a  Maxwellian  energy  spectrum,  because  in  this  case  the  ratios  scale 

exp  (-AE/kT) ,  where  AE  is  the  difference  in  excitation  energy  and  T  is 

the  (assumed  common)  temperature  at  which  the  lines  are  formed.  On  the 

other  hand,  photoionization  of  He  II  followed  by  recombination  to  excited 

states  with  subsequent  radiative  decay  to  ground  gives  line  ratios  in 

much  better  agreement  with  the  data  (Zirin) .  But  excitation  via  this  latter 

mechanism  implies  He  II  resonance  lines  are  formed  near  optical  depth  unity 

in  the  He  II  continuumand  hence  at  large  line-center  optical  depths.  If  this  were 

actually  the  case,  however,  the  emergent  line  profiles  should  exhibit 

a  central  self-reversal  (Milkey  1975),  and  they  do  not.  Instead  they  are 

observed  to  be  simple  Gausslans,  thereby  indicating  an  excitation  mechanism 

operative  at  small  line-center  optical  depths. 

The  above  dilemma  can  potentially  be  resolved  by  allowing  for 
collislonal  excitation  by  electrons  with  an  energy  spectrum  which  is  softer 
than  exponential,  such  as  those  calculated  here.  It  is- hoped  to  pursue 
this  and  other  aspects  of  the  helium  problem  in  a  later  publication. 
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V  .  Pi scussion 

a)  Implications  for  Transition-Region  Physics 

There  are  several  outstanding  questions  concering  the  low  TR 
(T  <  3  x  10^  K) ,  other  than  those  related  to  the  helium  spectrum 
discussed  in  §  IV,  for  which  electron  non-LTE  effects  hold  implications. 
Among  these  are  the  following. 

I)  Energy  Balance 

Energy  balance  in  the  low  TR  Is  widely  considered  to 

be  problematic  (e.g.,  Athay  1981,  Alvarez  1980,  Pallavicinl  et  al. 

1981,  Chiuderi  and  Riani  1974).  Broadly  stated,  the  problem  is 

that  empirical  temperature  profiles  are  inconsistent  with  known 

energy-transfer  mechanisms,  while  temperature  profiles  derived 

from  energy-balance  arguments  fail  to  reproduce  the  observations. 

The  shortcomings  of  classical  conduction-radiation  models 

are  especially  interesting  because  the  empirically  inferred  conductive 

5  2 

flux  at  the  top  of  the  quiet-Sun  TR  of  approximately  5  x  10  ergs/cm  -sec 

is  sufficient  to  balance  radiative  losses  from  the  transition 

region  and  upper  chromosphere  (Athay  1981).  There  are  two  difficulties.  On 

the  one  hand,  semi-empirical  analysis  of  EUV  line  intensity  data 

(e.g.,  Athay  1966,  Dupree  1972,  Kopp  1972:  see  also  Withbroe  1977) 

suggests  that  transition-region  temperature  gradients  are  roughly 

consistent  with  constant  conductive  flux  for  10^  K  £  T<10^  K,  but 

decrease  sharply  immediately  below  10^  K.  Since  the  classical 

5/2  dT  dT 

heat  flux  is  proportional  to  T  ,  a  decrease  in  gj-  implies  a 

large  conductive  energy  input  to  this  material  near  10^  K. 

As  recognized  by  Kuperus  and  Athay  (1967)  and  others,  however,  plasma  of 
solar  composition  at  10^  K  and  transition-region  pressures  cannot 


radiate  sufficiently  rapidly  to  dispose  of  the  incoming  energy, 
so  that  mass  motion  should  result.  Moreover,  as  a  result 
of  this  energy  deposition  near  10^  K,  there  is  an  insufficient 
conductive  energy  flux  into  the  upper  chromosphere,  where  the  bulk  of 
the  radiative  loss  actually  occurs.  Thus,  in  order  to  reconcile 
the  empirical  temperature  profile  with  energy  balance  one  requires 
an  energy  sink  near  10"*  K  and  an  energy  source  at  lower  temperatures 
(Athay  1981) . 

On  the  other  hand,  although  they  permit  conductive  energy 

in  amount  sufficient  to  balance  radiative  losses  to  reach  the 

upper  chromosphere,  (classical)  conduction-radiation  models  also 

prove  incompatible  with  energy  balance.  This  is  because  at  upper 

chromospheric  temperatures  the  temperature  gradient  required  to 

4 

carry  the  conductive  flux  is  so  large  (e.g.,  at  3  x  10  K, 

=  0.1  Km  if  9con<j  *  5  x  105  ergs/cm^-sec)  that  there  is 
insufficient  material  in  the  required  temperature  range  to  produce 
the  observed  radiation.  This  effect  was  encountered  in  §  III  of 
this  paper,  when  establishing  the  reference  temperature  scale  T  (z). 


The  above  difficulties  have  led  to  consideration  of 
alternatives  to  one-dimensional,  static,  conduction-radiation  models. 
However,  neither  alteration  of  the  geometry  (Gabriel  1976),  nor 
consideration  of  enthalpy  transport  by  steady  flows  (Pneuman  and 
Kopp  1978,  Chiuderi  and  Riani  1974)  have  significantly  changed  the 
situation;  to  account  for  the  observations  one  still  requires  an 
energy  sink  near  10^  K  and  an  energy  source  at  lower  temperatures. 
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Nov  consider  Che  possibilities  suggested  by  results 
reported  here  and  by  more  recent  calculations  (Paper  II).  First, 
it  is  plausible  that  rise  in  the  empirical  differential  emission 
measure  near  10^  K  is  due  to  electron  non-LTE  effects  on  ionization 
equilibria  and  line  emission,  as  illustrated  in  Fig.  7  for  magnesium 
and  discussed  in  detail  in  §  IV.  Lines  which  in  (electron)  LTE  are 
formed  over  a  relatively  narrow  temperature  range  near  10^  K  are 
perhaps  actually  formed  over  a  significantly  broader,  lower  temperature, 
and  hence  higher  density  range,  with  correspondingly  enhanced 
emission.  If  so,  there  need  not  be  a  sharp  decline  in  the  temperature 
gradient  immediately  below  10  K,  thus  obviating  the  need  for  an  energy  sink 
in  this  region.  This  view  is  supported  by  the  observational  evidence  dis¬ 
cussed  below.  Second,  based  on  recent  calculations  using  the  Bhatnagar, 

Gross  and  Krook  equation  (Paper  II),  I  suggest  that  the  conductive 
energy  flux  into  the  upper  chromosphere  is  carried  by  suprathermal  electrons 
with  velocities  extending  to  ID-30  times  local  thermal  velocities.  Since 
energy  deposition  by  suprathermal  electrons  is  due  primarily  to 
Coulour.b  collisions  with  background  electrons,  energy  deposition 
per  unit  pathlength  is  largely  independent  of  the  local  temperature 
structure.  Conductive  energy  transfer  via  this  mechanism  might 
therefore  overcome  the  small  scalelength  problems  inherent  to 
classical  conduction. 

To  summarize,  I  suggest  that  it  may  be  possible  to 
construct  kinetic  conduction-radiation  models  of  the  TR  which  are 
consistent  with  the  observations  and  with  energy  balance.  To  do 
so  requires  that  conduction  and  radiation  be  calculated,  and  observa¬ 
tions  interpreted,  with  proper  allowance  for  electron  non-LTE  effects. 


Within  the  context  "  this  paper  this  is,  of  course,  only  a  conjecture; 
detailed  calculations  are  required  in  order  to  reach  firm  conclusions. 
Moreover,  if  this  is  to  be  the  case,  the  role  of  energy 
transport  by  mass  motion  remains  to  be  clarified  (Athay  and  Holzer  1982). 
Nonetheless,  the  preliminary  results  reported  here,  the  simplicity 
and  economy  of  the  underlying  ideas  and  the  observational  evidence 
discussed  below  lend  some  weight  to  the  stated  conjectures, 
ii)  The  Schmahl  and  Orrall  (1979)  Observation 

The  above  authors  have  reported  evidence  for  strong 
continuum  absorption  of  EUV  line  photons  with  wavelengths  shortward 
of  912  &  everywhere  on  the  Sun's  disk.  The  lines  involved  had 
4.8  <  log  T  <5.6,  where  T  is  the  temperature  at  which  the 
standard  line  contribution  function  is  maximum.  As  discussed  by 
Schmahl  and  Orrall,  two  possible  explanations  are  i)  that  due  to  the 
inhomogeneous  nature  of  the  atmosphere  there  is  a  large  amount  of 
cool,  neutral-hydrogen-containing  material  overlying  the  EUV -emitting 
regions  (most  likely  in  their  opinion),  or  ii)  that  there  may  be  a 

significant  amount  of  neutral  hydrogen  (and  helium)  distributed  more 
or  less  uniformly  with  the  EUV-emitting  ions. 

I  note  that  provided  the  additional  effects  discussed  in 
§  IV d  do  not  qualitatively  change  the  ionization  balance,  the  second 
of  the  above  alternatives  is  what  one  would  expect  on  the  basis  of  the  ideas 
discussed  in  this  paper.  Ionic  stages  which  in  electron  LTE  exist  only  at 

relatively  high  temperatures  (2-3  x  1C5  K) ,  may  actually  have  significant 


populations  throughout  the  upper  chromosphere  where  there  is  a  non- 
negligible  concentration  of  neutral  hydrogen. 
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What  speaks  against  this  Intepretation  is  that  the  absorption  appears 
to  be  as  pronounced  in  cell  interiors  as  in  network  regions.  However,  the 
magnetic  field  configuration  in  both  network  and  cell  interior  regions 
is  perhaps  too  uncertain  for  this  objection  to  be  decisive. 

iii)  Temperature  Plateaus 

Empirical,  one-dimensional  models  of  the  upper  chromosphere 
require  substantial  amounts  of  material  at  temperatures  near  20,000  K 
in  order  to  account  for  the  observed  resonance  line  and  continuum 
spectrum  of  hydrogen  (Vernazza  e£  al.  1976,  1980,  Basri  et  al .  1979, 
Gouttebroze  et  al^.  1979).  Although  the  existence  of  such  plateaus  can 
plausibly  be  attributed  to  inhomogeneities  in  the  solar  atmosphere,  within 
the  context  of  one-dimensional  models  they  are  problematic  in  two 
respects  (Avrett  1981).  First,  it  is  difficult,  from  an  energy-balance  stand¬ 
point,  to  identify  processes  responsible  for  their  formation,  and 
second,  spectrum  synthesLs  for  elements  other  than  hydrogen  often 
require  plateaus  with  disparate  properties.  For  example,  Lites  et  al . 

(1978)  find  that  their  best  fit  to  center-to-limb  measurements  of  C  II 
resonance  lines  at  133.5  nm,  as  well  as  hydrogen  Ly-a  and  Lyman 
continuum  intensities  requires  a  plateau  at  16,500  K  with  about  25% 
more  material  than  the  one  at  20,000  K  required  by  Vernazza  et  al. 

(1976)  based  on  hydrogen  data  alone  (see  Vernazza  et^  al^.  1981). 

Loulergue  and  Nussbaumer  (1976)  note  that  the  observed 

1977  and  X1176  emission  from  C  III  can  be  explained  by  a  plateau  at 

10  -3 

30,000  K  with  an  electron  density  n  ■  10  cm  and  a  thickness  of  800  km 
along  the  line  of  sight.  This  is  especially  noteworthy  in  the  present 
context  because  standard  calculations  (Jordan  196y)  indicate  that  C  III 
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has  its  maximum  abundance  near  60,000  K.  Finally,  Feldman,  Doschek  and 
Patterson  (1976)  find  that  X 1207  emission  from  Si  III  is  compatible  with 
a  plateau  at  36,000  K  with  thickness  840  km  in  cell  boundary  regions  and 
380  km  in  cell  interiors. 

I  suggest  that  the  plateaus  just  mentioned  are  perhaps  artifacts 
of  the  LTE  assumption;  i.e.,  they  are  required  in  present  empirical 
models  in  order  to  compensate  for  underestimating  collisional  excitation 
and  ionization  rates. 

iv)  Non-local  Thermodynamics 

It  follows  from  results  shown  in  5  IV  that  the  detailed 
shape  of  the  electron  distribution  function  in  the  low  TR  is  dependent  on 
the  global  temperature  and  density  structure  of  the  overlying  atmosphere. 

An  implication  of  this  is  that  the  ionization  equilibrium  of  the  elements, 
optically  thin  radiative  losses  and  energy  transport  by  thermal  conduction 
acquire  a  dependence  on  the  global  structure  of  the  upper  TR  and  corona. 

This  siutation  is  in  sharp  contrast  to  that  which  attains  in  electron-LTE, 
where  each  of  the  above  processes  is  a  function  of  the*local  thermodynamic 
variables  and  perhaps  their  gradients.  Clearly,  breakdown  of  the  local 
Maxwellian  approximation  in  the  solar  TR — and  by  implication  in  the  transition 
regions  of  other  stars — leads  to  a  significant  increase  in  the  level 
of  complexity  required  for  successful  spectroscopic  diagnostic  work  or 
theoretical  modeling. 

If  future  calculations  give  fractional  abundance  curves  similar 
to  the  preliminary  results  for  magnesium  shown  in  Fig.  7,  two  added 
complications  which  may  then  arise  are  the  necessity  of  considering  charge 
transfer  r ^combination  reactions  in  the  ionization  balance  and  finite 
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optical-depth  effects  in  resonance  lines  which  are  optically  thin  according 
to  standard  ionization  equilibrium  calculations. 

Finally,  I  note  that  because  the  calculated 
electron  energy  spectrum  has  a  relatively 

weak  temperature  (spatial)  dependence  in  the  low  TR  (Fig.  3),  ionization 
nonequilibrium  effects  due  to  flows  through  this  region  (Joselyn  et  al. 

1978,  Raymond  and  Dupree  1978)  may  be  less  significant  than  presently 
thought . 

v)  Other  Effects 

The  most  straightforward  test  of  the  ideas  discussed  here  lies 
in  comparison  of  predicted  and  observed  EUV  emission-line  profiles  and 
intensities.  However,  there  are  two  additional  effects  associated  with 
a  nonthermal  electron  spectrum  which  may  be  observable. 

As  mentioned  in  §  II,  it  is  well  known  that  electron  distributions 
of  the  sort  calculated  here  imply  an  enhanced  level  of  electron  density 
fluctuations  at  frequencies  near  the  local  plasma  frequency.  (If  the 
electron  temperature  exceeds  the  ion  temperature,  ion  density  fluctuations 
will  also  be  excited.)  In  addition  to  their  ability  to  scatter  fast 
electrons  (Tidman  and  Eviatar  1965) — an  effect  which  should  be  included 
in  future  calculations — Tidman  and  Dupree  (1965)  and  others  have  shown  that 
such  density  fluctuations  also  lead  to  a  significant  collective  contribution 
to  the  bremsstrahlung  emission  at  the  local  plasma  frequency  and  its  first 
harmonic.  For  the  TR  the  relevant  frequency  range  is  roughly  .1  to  10  GHz. 
t  find  that  for  a  model  atmosphere  consisting  in  an  isobaric,  constant 
conductive  flux  TR  (10^  <_  T  <_  2  x  10^  K,  NgTe  *  6  x  10^  cm  q  *  6  x  10^ 
cm  “  sec  ),  together  with  an  isothermal,  hydrostatic  corona 


ergs 


(Tcor  *  2  x  10^  K) ,  the  optical  depth  is  less  than  unity  for  frequencies 
v  >  7  GHz  and  .3  GHz  for  0;^  and  2-a)p  emission,  respectively.  Although  . 

these  numbers  are  quite  model  dependent,  they  indicace  that  it  may  be 

possible  to  observe  enhanced  (o  and  2-w  radiation  from  the  TR. 

P  P 

Henoux  £t  al.  (1982)  have  observed  linear  polarization  in  the 
chromospheric  S  I  1437A  line  during  the  soft  x-ray  phase  of  a  flare.  The 
measured  direction  of  polarization  suggests  that  collisional  excitation 
by  an  anisotropic  electron  distribution,  such  as  those  expected  in  a 
plasma  carrying  a  heat  flux,  is  the  mechanism  responsible  for  producing 
the  polarization.  If  higher  temperature  lines  suitable  for  polarization 
measurements  could  be  identified,  and  measurements  made,  they  might  provide 
a  valuable  constraint  on  future  calculations  of  the  electron  distribution 
function. 

b )  Previou s  Work 

Spicer  (1978)  has  noted  that  the  large  heat  flux  inferred 
from  empirical  models  of  the  solar  TR  imply  a  return  current  electric 
field  which  is  of  the  order  of  the  Dreicer  field.  H°  suggests  that  as 
a  result  there  occurs  a  "boiling  off"  of  electrons  in  the  tail  of  the 
electron  distribution  and  that  these  electrons  represent  a  non-negligible 
fraction  of  the  total  energy  input  to  the  corona. 

This  idea  is  incorrect.  Present  results  show  that  in  a 
temperature  gradient  fast  electrons  diffuse  down  the  gradient  against 
the  action  of  the  electric  field.  Spicer's  error  lies  in 
accepting  conventional  ideas  concerning  electron  runaway  in  a  homogeneous 
plasma  immersed  in  a  dc  electric  field  (Dreicer  1959).  Reference  to  Eq.  41 
shows  that, in  an  inhomogeneous  plasma,  gradient-related  terms  dominate 
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electric  field  terms  by  a  factor  of  (— )2  in  the  kinetic  equation.  The 

Vth 

latter  are,  therefore,  unimportant  at  suprathermal  velocities. 

Pursuing  ideas  (Shoub  1976)  made  available  to  him,  Roussel -Dupre  (1980 
has  examined  the  effect  of  an  electron  distribution  function  possessing  a 
nonthermal  tail  on  heat  conduction  and  on  the  ionization  equilibrium  of  several 
elements,  because  there  is  no  physical  basis  for  his  choice  of  distribution 
function,  however,  meaningful  conclusions  cannot  be  drawn  from  his  analysis. 

Gurevich  and  Istomin  (1979)  have  attempted  a  perturbation  analysis 
of  the  linearized,  high  velocity  form  of  the  Fokker-Planck  equation,  treating 
a  -  as  a  small  parameter.  Because  they  resort  to  quite  severe 

idealizations  (e.g.,  infinitely  steep  temperature  gradients), 
their  results  appear  to  be  only  of  qualitative 
significance.  Nevertheless,  they  do  point  out  classical  heat  conduction 
results  are  likely  to  be  invalid  for  quite  snail  values  of  a. 

c)  Summary 

In  the  foregoing  I  have  shown  that  the  fundamental  assumption  that  in  a 
weakly  inhomogeneous  plasma  free  electrons  maintain  nearly  a  local-Maxwellian 
velocity  distribution  function  is  invalid  throughout  the  solar  transition 
region  and  upper  chromosphere.  I  have  examined  several  implications  of 
the  breakdown  of  this  assumption  and  have  suggested  that  a  number  of  out¬ 
standing  puzzles  related  to  the  low  transition  region — including  energy 
balance,  the  anomalous  resonance  line  spectrum  of  He  I  and  He  II,  the 
anomalous  continuum  absorption  by  neutral  hydrogen  and  the  seeming  need 
for  plateaus  in  empirical  one-dimensional  temperature  profiles  of  the 
upper  chromosphere — are  potentially  resolvable  by  abandoning  the  electron-LTE 
hypothesis.  In  addition,  I  have  identified  several  phenomena— non-local 


heat  conduction,  enhanced  electron~ion  excitation,  ionization  and  - 
dielectronic  recombination  rates,  charge  transfer  recombination  on 
neutral  hydrogen  and  helium  and  possible  plasma  collective  effects — 
which  may  play  a  role  in  future  development  of  this  subject.  These 
advances  were  made  possible  by  development  of  an  effective  numerical 
algorithm  for  solving  kinetic  equations  of  the  Fokker-Planck  type. 
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Table  1. 

a 

Empirically  Derived,  Homogeneous  Quiet-Sun  Transition  Region  Model 


z 

Km 

log  To 

°K 

no 

-3 

cm 

dT 

_ 0 

dz 

°K/Km 

A 

Km 

a 

_1 

5c“a4 

2006 

4.8 

9.5(9) 

4.3(3) 

3. 3(-2) 

2. 2 (-3) 

4.6 

2011 

5.0 

6.0(9) 

1.1(4) 

1. 2(-l) 

1.3 (-2) 

2.9 

2016 

5.2 

3.8(9) 

1.1(4) 

4.6 (-1) 

3.2(-2) 

2.4 

2031 

5.4 

2.4(9) 

3.6(3) 

1.7 

2.4 (—2 ) 

2.5 

2512 

5.8 

9.5(8) 

1.1(2) 

26.4 

4.6 (— 3 ) 

3.8 

a 


b 


A.  K.  Dupree  (1972) 

NOTE  -  Numbers  in  parenthesis  denote  multiplicative  powers  of  10 


Table  2a. 


First  Two  Angular  Moments  of  the  Distribution  Function  for  Three  Coefficient 

Iterations  for  the  Case  P  *  P  at  a  Location  Where 

o 

t,  =  27,  T  -  19.900K,  a  -  3.8  x  lO-2 


♦  (5.  t,)/**U) 


[teration 

Vno. 
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1  2 

«)“ 
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1.00  .992 
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3.73 

(-2) 

-2.36 

(-2) 

3.02 

(-2) 

-3.15 

(-2) 

-3.01 

(-2) 

1.13 

-1.18 

-1.16 

8.94 

(1) 

-9.43 

(1) 

-9.25 

(1) 

1.21 

(5) 

-1.27 

(5) 

-1.25 

(5) 

C  =  1/2  J  £,  t.)  dy;  $  =  3/2  I  y  5,  t.)  dy 

-1  -1 


First  iterate  corresponds  to  MFPA. 
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Table  2b. 

Calculated  Dimensionless  Heat  Flux  for  Three  Coefficient  Iterations 

At  Same  Location  as  in  Table  2a 


Iteration3 

q°*C 

1 

—  1 . 04  (— 1 ) 

2 

-9. 46 (-2) 

3 

-7.47(-2) 

First  iterate  corresponds  to  MFPA. 


b  q  =  /  u  du  J  C5  ♦  <u.5.t1) 

-1  0 

c  Spitzer-Harm  value  for  q  is  -8.06 (-2) 
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Table  4. 

Kinetic  Theory  Values  of  Reduced  Thermodynamic  Variables 


ft 

Tq(°K) 

n/n 

0 

T/T 

0 

P  -  P 

0 

A  i 

N 

V 

% 

H 

1.12(4) 

1.011 

1.008 

-0.002 

1.016 

3.15(4) 

1.006 

0.999 

-0.002 

1.006 

6.27(4) 

1.004 

0.999 

-0.002 

1.003 

1.17(5) 

1.003 

0.999 

-0.001 

1.003 

2.39(5) 

1.002 

0.998 

-0.001 

1.003 

% 

6.93(5) 

1.001 

0.999 

-0.001 

1.001 

,  P"/ 
<  £  >  i 


1.016  1.12(4)  1.020  1.026  0.004  1.027 

1.006  3.23(4)  1.012  1.007  -0.002  1.015 

1.003  6.55(4)  1.009  1.002  -0.003  1.010 

1.003  1.10(5)  1.007  1.000  -0.002  1.007 

1.003  2.35(5)  1.005  0.999  -0.002  1.004 

1.001  6.90(5)  1.003  0.998  -0.001  1.002 


Table  5 


;  ; 

3 


Normalized  Isotropic  Component  of  Distribution  Function 

4>°(S,To)  /**U) 


5 

T  -  1.5 

X  io4k 

T  -  3.1 

X  104K 

T  -  1.57 

X  105K 

T  -  6. 

92  X  105K 

o 

0 

0 

o 

P 

P  n 

P 

P  /2 

P 

P  12 

P 

P  /2 

o 

O 

0 

0 

0 

0 

o 

0 

0.00 

0.99 

0.98 

0.99 

0.98 

1.00 

1.00 

1.00 

1.00 

0.28 

0.99 

0.98 

1.00 

0.98 

1.00 

1.00 

1.00 

1.00 

0.57 

0.99 

0.99 

1.00 

0.99 

1.00 

1.00 

1.00 

1.00 

0.85 

1.00 

1.00 

1.00 

1.00 

1.00 

1.13 

1.00 

1.00 

1.00 

1.00 

1.00 

1.41 

1.00 

1.00 

1.00 

1.00 

1.00 

0.99 

0.98 

1.00 

0.98 

1.00 

1.00 

1.00 

1.98 

0.98 

0.94 

0.99 

0.94 

1.00 

0.98 

1.00 

1.00 

2.26 

0.97 

0.92 

0.98 

0.92 

1.00 

0.97 

1.00 

0.99 

2.54 

1.02 

1.08 

1.00 

1.08 

1.00 

0.99 

1.00 

0.99 

2.82 

1.3 

1.8 

1.2 

1.8 

1.03 

1.2 

1.01 

1.03 

3.  10 

2.4 

4.1 

1.8 

4.1 

1.21 

1.8  - 

1.04 

1.2 

3.39 

6.1 

1.5(1) 

4.0 

1.3(1) 

1.8 

3.9 

1.2 

1.8 

3.67 

2.1(1) 

6.1(1) 

1.2(1) 

3.4(1) 

3.74 

1.2(1) 

1.6 

3.9 

3.95 

9.3(1) 

3.2(2) 

5.0(1) 

1.6(2) 

1.1(1) 

4.9(1) 

3.1 

1.1(1) 

4.24 

5.2(2) 

2.0(3) 

2.6(1) 

9.5(2) 

4.5(1) 

2.5(2) 

8.0 

4.4(1) 

4.52 

3.6(3) 

1.6(4) 

1.7(3) 

6.9(3) 

2.4(2) 

1.6(3) 

2.9(1) 

2.2(2) 

4.80 

3.0(4) 

1.5(5) 

1.4(4) 

6.1(4) 

1.6(3) 

1.3(4) 

1.4(2) 

1.3(3) 

5.08 

3.1(5) 

1.6(6) 

1.3(5) 

6.5(5) 

1.4(4) 

1.3(5) 

8.8(2) 

9.8(3) 

5.36 

3.8(6) 

2.2(7) 

1.6(6) 

8.3(6) 

1.5(5) 

1.6(6) 

6.8(3) 

8.2(4) 

5.65 

5.7(7) 

3.5(8) 

2.3(7) 

1.3(8) 

1.9(6) 

2.2(7) 

6.3(4) 

7.8(5) 

5.93 

1.0(9) 

4.6(9) 

3.9(8) 

2.3(9) 

3.0(7) 

3.8(8) 

6.3(5) 

8.2(6) 

1 


Figure  Captions 
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Figure  L.  A  plot  of  the  pitch  angle  dependence  of  the  reduced  distribution 

function  at  a  fixed  velocity  and  location,  illustrating  the  conver¬ 
gence  properties  of  the  iterative  scheme  discussed  in  the  text. 

Curves  arc  labeled  by  iteration  number.  Thus,  curve  "0"  is  the  initial 
estimate,  curve  "1"  is  the  result  of  the  first  forward  (y>0)  integra¬ 
tion,  curve  "2"  is  the  result  of  the  first  backward  (y«0)  integration, 
ate. 


Figure  2.  Isotropic  part  of  the  reduced  distribution  function  versus  dimension- 
le 

r2 


less  energy  at  several  heights  for  P  =*  P^.  Abscissa; 


2M  (Zy  Ordinate;  4>°  (5,  TQ)  -  1/2  ^  (y,  5,  TQ)  dy.  Solid 


lines  are  calculated  results.  Dashed  line  is  a  local  Maxwellian. 
Solid  curves  are  labeled  with  the  value  of  at  the  corresponding 


height . 


Figure  3.  Isotropic  part  of  the  dimensional  distribution  function  versus  velocity 

g 

measured  in  units  of  10  cm/sec  for  the  case  P  ■  Po/2’  Das^e^  lines  are 
local-Maxwellian  distributions;  solid  lines  are  calculated  results. 


Figure  One-dimensional  reduced  distribution  function  of  electron  velocities 
along  the  gradient  direction  at  several  heights,  including  upper  and 
lower  boundaries.  is  positive  in  the  direction  of  increasing 

temperature.  Solid  lines  represent  calculated  results.  Dashed  lines 
are  local  Maxwellians. 


Figure  5.  Normalized  pitch  angle  dependence  of  the  distribution  function  for  a 
sequence  of  velocities  at  a  location  near  the  bottom  of  the  slab 
(IA  *  17.500K).  u  is  the  positive  in  the  direction  of  increasing  temper¬ 


ature.  Note  ordinate  scale  changes  between  panels. 


Figure  6.  Contour  plot  of  the  reduced  distribution  function  for  P  •  Pq  and  at 
the  same  location  as  in  Fig.  5.  I,-  •?*  '-ha  component  of  the  dimen¬ 
sionless  velocity  alone  the  gradient  direction  and  is  positive  in  the 
direction  ol  increasing  temperature. 

Figure  7.  Normalized  pitch  angle  dependence  of  the  distribution  function  at  the 
lower  boundary.  The  boundary  condition  is  4*  (y,  €»  T  *  0)  »  y*  (£) 
for  y>0. 

Figure  8.  Fractional  relative  abundance  of  the  first  six  ionization  stages  of 

magnesium  versus  temperature  for  P  *  P^  and  Fq^2.  Solid  lines  are 


for  calculated  collision  rates.  Dashed  lines  are  for  Maxwellian  rates. 


Figure  5 


Figure 
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Appendix 

Here  I  give  heuristic  arguments  supporting  the  step-size  estimates  (49a, b). 

Consider  first  the  case  of  explicit  differencing.  In  their  book  Richtmyer 
and  Morton  (1968,  p.  207)  use  a  standard  Von  Neumann  stability  analysis  to  show, 
under  the  condition  that  the  coefficients  be  slowly  varying,  that  if  Eq.  (47) 
is  differenced  using  a  forward  (explicit)  time  difference  and  centered,  second- 
order  space  differences,  the  stability  condition  is 


At 


( 


+ 


(Al) 


This  is  a  local  condition  in  that  it  must  be  satisfied  for  all  x,  y  in  the 
region  of  interest. 

Disregarding  the  fact  that  the  coefficients  in  Eq.  45  are  not  slowly 
varying,  the  analogous  condition  for  stability  of  an  explicit  differencing  scheme 
for  that  equation  (using  centered  second-order  u  and  £-dif ferences)  is 


-  /  A  C  \  1 

"  (co2  u,)2r2 


(A2) 


2  -3 

From  Eq.  (43)  we  see  th  t  in  the  MFPA,  A  -*•  and  C  -*■  Z  as  ;  0.  Since 

j 

the  smallest  values  of  u  and  5  on  the  mesh  are  —  and  AC,  (A2)  implies 


A-  <  0  j(.p)  (.-.02  ^ 


(A3) 


which  is  (49a)  of  the  text.  Due  to  the  rapid  variation  of  the  coefficients,  (A3) 
is  most  probably  incorrect  as  a  rigorous  stability  criterion.  Nevertheless, 


it  is  probably  not  very  wrong  and  thus  indicates  that  the  singular  nature  of  the 
coefficients  in  the  Fokker-Planck  equation  near  5*0  leads  to  excessively  small 
step  sized  restrictions  for  explicit  spatial  differencing  schemes. 


Now  consider  application  of  ADI-dif ferencing  (Richtmyer  and  Morton,  p.  211) 
to  Eq.  (48b)  in  the  MFPA,  so  that  by  (39b)  the  mixed-derivative  term  is  absent. 
In  this  case  (48b)  can  be  written 


If  -  %  -  i+(,)  +  i  « 


where  for  5  t  0,  5  ,  and  L  are  difference-matrix  representations  of  the 

max’  «5  2  r 

operators  yr-  ja*3^  +  £d*  +  V  (^| - 8)  J  3^  +  P*  +  y  auC  ^  and 


operators  yr  |a*3^  +  £d*  +  y  (5i| - g)J 

|C”  3UU  +  [e*  +  B]  respe, 


respectively.  The  vector  q  denotes  the 


+  -  ' 

term  S  i  in  (48b)  and  is  considered  known.  Elements  of  L  and  L.  for  5  *  0,5 
»  —  «y  »5  n 

correspond  to  the  difference  form  of  the  side  conditions  (44)  and  (42)  respec¬ 
tively.  Using  the  notation  »  4>(i£.x)  and  omitting  the  +  superscript  from  b+, 
the  ADI  scheme  for  (A4)  may  be  written 

c  -  c  T  ,1+1/2  ,  T  ,i  ,  i 

- - : — )r.  *  *  L.  0  +  I.  <(>  +  q 

_t/2  *5  —  *y  1  2 


il+1  -  i1+1/2  _  ,1+1/2  ,1+1  1+1 

- - 1772 - h  i  +  K  1  +  3. 


i  1  /  2 

Upon  eliminating  61  I  find  that 


- rr-*—  -  4-  (L,  +  l  )  (  ^  +  *i+1  \  +4  Iq1  +  QX+l\ 


At 
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I  observe  that  the  first  two  bracketed  terms  on  the  right  side  of  (A6)  are 
the  result  for  a  fully  implicit  (hereafter  FI)  step,  centered  at  i  +  1/2.  The 


last  term  on  the  right  in  (A6)  is  thus 


3d  3© 


In  usual  applications 


FI 


of  the  ADI  method  (slowly  varying  coefficients)  this  difference  is  of  order 
2 

(At)  ,  and  is  acceptably  small  because  the  truncation  error  of  the  FI  step  is 

2 

itself  of  order  (At)  .  In  our  case,  however,  this  term  is  large  at  small 
velocities. 

To  see  this  we  make  use  of  the  fact  that  for  conditions  of  interest  the 
computed  solution  is  close  to  the  (1953)  Spitzer-HMrm  solution  at  small 
velocities;  i.e., 

©  (p.,5,t)  s  <p  (?)  +  2u  a  (t)  d1  (C) 

where 

d1  (5)  -  |dt  (5)  -  0.35  D£  (?)  |  d* 

and  D  and  D  are  the  functions  tabulated  by  Spitzer  &  Harm.  I  find  that 

1  Cj 

1 

:  (?)  »  0  (?  )  as  ?  +  0  where  0<o << 1 .  Using  the  definitions  given  above 

for  the  operators  L  and  L-  I  now  find  that 


L„  L 


(-■l+l  -  -  o 

da  1 

\  At  )  5-0 

dT  3-6 

L  u?  J 

(A7) 


3  © 
3t 


Since  4^  =  0(a3),  the  requirement  that 

dT  i “ 1  ADI  W‘FI 

with  help  of  (A6)  and  (A7),  to  the  condition  that 


3© 

3t, 


<<  ©  for  all  u,C  leads. 


c»At  <  o  £a  u  (a-;)3J 


1/2 


which  is  (49b)  of  the  text. 
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